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ABSTRACT 


Earth imaging satellites have typically been large systems with highly 
accurate and expensive sensors. With the recent push for Operationally 
Responsive Space, Earth imagining has become potentially achievable with 
small and relatively inexpensive satellites. This has led to the research currently 
underway to develop very small, low-cost imaging satellites that can produce 
useful Operational-level and Tactical-level imagery products. This thesis 
contributes to that effort by developing a detailed design for the attitude 
determination system for a tactically useful earth-imaging nano-satellite. Tactical 
Imaging Nano-sat Yielding Small-Cost Operations and Persistent Earth-coverage 
(TINYSCOPE) is an ongoing investigation at NPS, concerning using a nano¬ 
satellite, based on the CubeSat standard, to achieve Earth imaging from LEO 
orbit. 

A detailed design of the attitude determination system includes sensor 
selection and characterization, as well as high fidelity simulation via 
MATLAB®/Simulink®. The attitude determination system is based on an 
Extended Kalman Filter using multiple sensor types and data rates. The sensors 
include a star tracker, sun sensor, gyroscope, and magnetometer. 
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I. BACKGROUND 


A. BRIEF HISTORY 

Since the very early days of satellite imagery, resolution and pointing 
accuracy have been key parameters of the satellite design. Corona, America’s 
first Earth Imaging satellite [1], had a resolution of 8-1 Om. This was eventually 
improved to 2-4m resolution through improved cameras and lower altitudes. As 
subsequent generations of military Earth imaging satellites were being 
developed, the resolution continued to improve, which in turn lead to the need for 
higher accuracy pointing. The new satellites used ever-increasingly complex and 
expensive attitude determination systems (ADS) to provide the ability to meet 
these new requirements. 

Today, commercial imagery satellites are going through the same trend in 
resolution. They have improved in the past decade, from relatively low resolution 
at about 5m to about 0.6m [2]. At the same time, these commercial satellites are 
getting smaller and keeping costs to a minimum. These two facts present new 
challenges to the ADS. The traditional sensors are too big and expensive to be 
practical for the new satellites. New sensors have been shrinking due to smaller 
electronics; however, this miniaturization is not keeping pace with the rapidly 
increasing demands on the ADS. 

The latest Earth imagers take the trend of shrinking satellites to a new 
level. They are pico-satellites based on a new standard call CubeSats. They 
push the bounds of small size and low cost. Developed to provide cheap access 
to space and used primarily by universities, these satellites are measured in the 
tens of centimeters, weigh between one and ten kilograms, and usually cost in 
the tens of thousands of dollars. CubeSats may become a viable Earth imaging 
platform to complement large and very high-resolution existing assets. 
Therefore, the ADS is becoming a much more important part of the design. High 
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accuracy attitude determination is still in its infancy for CubeSats. This thesis will 
attempt to develop an ADS that meets the needs of this new class of satellite. 

B. CUBESAT STANDARD 

In 1999, California Polytechnic University (Cal Poly), San Luis Obispo, and 
the Space Systems Development Laboratory (SSDL) at Stanford University 
collaborated to develop a standard for a pico-satellite design [3]. The goal was to 
enable rapid development and launch of satellites in the most cost effective 
manner. The target audience was universities and research institutes and there 
are currently over 100 universities, high schools, private firms, and government 
agencies developing CubeSats. The concept would enable groups to test new 
and innovative hardware and software on actual satellites instead of relying on 
simulations. The standardized size and mass would enable a standardized 
launcher to be developed that would be suitable for secondary or tertiary launch 
opportunities. It also enables the development of large numbers of experienced 
Astronautical Engineers. A student can now, in principle, with sufficient available 
funding and professional support, design, build, test, and launch a satellite during 
his undergraduate or graduate school period. 

The CubeSat standard also allows for expansion from a single cube (1U). 
Multiple cubes can be attached together to form larger CubeSats like the 3U (3 x 
1 cubes). The number of cubes attached together is really only limited by the 
launching mechanism. Launchers that can handle 5U (5 x 1 cubes) and 6U (2 x 
3 cubes) are under development. These larger CubeSats can accommodate 
payloads that are more complex. The flexibility of this standard is also reflected 
in the fact that it is an open standard. Not just the original creators, but also the 
entire community of users continuously review the standard. An open standard 
can evolve to meet the needs of the community. The latest revision to the 
CubeSat mechanical requirements (Rev 12) are summarized in Table 1 and 
Figure 1. Another important design consideration for flexibility and utility was the 
electronics that would make up the working parts of the CubeSat. The sizing of 
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the CubeSats allows for commercial electronics in the PC/104 form factor and 
commercial solar cells to be used. This opens a large array of previously existing 
systems to be used in the CubeSats. On the same point, it allows for easy 
development of new electronic boards and instruments using this form factor. 



Figure 1. 1U CubeSat Side Definition. From [4]. 


# 

Requirement 

Value Unit 

1 

X and Y dimensions 

IOOiO.1 mm 

2 

Z dimensions (1U) 

113.5±0.3 mm 

3 

Z dimensions (3U) 

340.5±0.3 mm 

4 

Maximum component protrusion from X or Y side 

6.5 mm 

5 

Mass (1U) 

1.33 kg 

6 

Mass (3U) 

4.0 kg 

7 

Center of gravity and geometric center difference 

20 mm 


Table 1. CubeSat Mechanical Specifications. After [4]. 
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C. CUBESAT DEPLOYMENT 

Cal Poly has also created a deployment standard Poly Picosatellite Orbital 
Deployer (P-POD). This is a box structure that houses the CubeSat during 
launch (see Figure 2). It provides a standard attachment point to the launch 
vehicle and a deployment mechanism for the CubeSat. It protects the CubeSat, 
launch vehicle and the primary payload. This is particularly important because 
the owners of the primary payload can be very conservative and want 
guarantees that their very expensive satellite is not damaged by the secondary 
payload. To ensure this, the P-POD has been tested to very high standards and 
proven through numerous launches that it does this job very well. It has also 
proven to be very flexible. It is compatible with a number of launch vehicles (see 
Table 2), and any CubeSat can be launched from it. This means that if the 
manifested CubeSat cannot launch for some reason, another one can replace it 
very quickly and easily. This is because any CubeSat made to the CubeSat 
Standard fits in a P-POD and the launch vehicle only cares about the P-POD. 
This combination of fitting many launch vehicles and flexible payload manifesting 
enables much easier access to space than ever before. Currently the P-POD 
can hold up to three 1U or one 3U CubeSats. There is currently research 
underway to develop a 5U P-POD and a 2 x 3 or 6U P-POD at Cal Poly. 


Delta 11“ 

Peeasus 

Delta r\" — SAM 

Taurus 

Delta TV - ESP A- 

Falcon I 

Atlas — ESPA*'" 

Rockot (Russian) 

Space Shuttle 

Dnepr’^ (Russian) 


Table 2. Launch vehicles compatible with CubeSat P-PODs. 

From [3]. 
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Figure 2. CubeSat P-POD Unit. From [3]. 

The Naval Postgraduate School (NPS) is currently working with Cal Poly 
to meet the need for higher capacity CubeSat launches. The collaboration is 
developing a high capacity CubeSat launcher that will be designed to attach to 
the Evolved Expendable Launch Vehicle Secondary Payload Adapter (ESPA). 
The NPS CubeSat Launcher (NPSCuL) will be a collection of 10 P-PODs (see 
Figure 3) with coordinating electronics for CubeSat deployment. NPSCuL will be 
able to deploy a combination of 1U, 3U, 5U, or 6U CubeSats. 

Cal Poly is not the only developer of CubeSat launchers; several other 
CubeSat launchers have been developed and launched [5]. Tokyo Institute of 
Technology’s Lab for Space Systems (LSS) has developed the Tokyo Pico- 
satellite Orbit Deployer (T-POD). Germany’s Astrofein has developed the Single 
Pico-Satellite Launcher (SPL). Both of these are 1U deployers. The most used 

5 


Cubesat launcher, other than the P-POD, is the experimental Push Out Deployer 
(X-POD) developed by University of Toronto Institute for Aerospace 
Studies/Space Flight Laboratory (UTIAS/SFL). At least six X-PODs have 
launched to date. They come in a variety of configurations including 1U, 3U, and 
the X-POD DUO holds a satellite 20 x 20 x 40 cm. NASA has even started 
developing a CubeSat launcher to help further develop their CubeSat program. 
They are specifically designing a 6U or “Six Pack” deployer. 



Figure 3. NPSCuL Model. From [6]. 

D. SURVEY OF CUBESAT ATTITUDE DETERMINATION SYSTEMS 

The majority of the CubeSat mission that have launched so far have some 
sort of attitude sensing (see Table 3). However, these sensors are usually low 
accuracy sensors like magnetometers or sun sensors. This is due to several 
factors. First, most of the missions flown so far do not have active attitude 
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control, which leads to the assumption that attitude determination is also of little 
use. Second, most of the CubeSats with active 3-axis control typically use only 
magnetorquers. These provide a simple low-cost control but are also low 
accuracy; again, no need for high accuracy attitude determination. Finally, the 
high accuracy ADS is much more complex and expensive to implement. A very 
high accuracy sensor (<100 arcsec), like a star tracker, has yet to be flown. 
There are several currently in development specifically for CubeSats, but only 
one currently commercially available that can fit in a CubeSat. 
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Tokyo bistitiite of Technoloav 

CUTE-I 

2003 

K' E MS Gyro. 2-Axis Accetercmeter. Sun Sensor 

Universitv of Tokyo 

XJ-IV 

2003 

None 

University of Toronto 

CanX-1 

2003 

Horizon sensor & Star tracker. Magnetometer. GPS 

Technical University of Denmark 

DTUsat 

2003 

MEMS sun sensor. Maiyielomeler 

Alborg University 

AAUCubesat 

2003 

Sun sensor. Magnetometer 

Stanford University and Quakesat U.C 

(XakeSat 

2003 

Magnetometer 

Nonwegian U of Science and Technology 

NCube2 

2005 

Magnetometer 

University of Wurzburg 

UWE-1 

2005 

N/A 

Universitv of Tokyo 

XJ-V 

2005 

Camera 

Tokyo Institute of Technology 

CUTE 1.7+ APD 

2006 

Gyro. Magnetometer. Sun sensor. Earth sensa 

Universitv of lllionis 

ION 

2006 

Sun sensors, magnetometer 

University of Arizona 

Sacred 

2006 

NIA 

University of Kansas 

KUTEsatPathfintk 

2006 

Magnetometer, sun sensors 

Cornell Universitv (New York state) 

ICC Cube 

2006 

Magnetometer. GPS receiver 

University of Arizona 

RINCON 1 

2006 

N/A 

Nihon University 

SEEDS 

2006 

Gyro sensor, and magnetometer 

Hankuk Aviabon Universitv 

HAUSAT 

2006 

Sun sensor, gps receiver 

Norwegian U of Science and Technology 

Ncube 1 

2006 

Magnetometer 

Montana State University 

MEROPE 

2006 

N/A 

Aerospace Corporation 

AeroCube-1 

2006 

N/A 

California Polytechnic Institute 

CP2 

2006 

mavielameter 

California Polytechnic Institute 

CPI 

2006 

Sun sensor 

Cornell Universitv (New York state) 

ICE Cube 2 

2006 

Magnetometer. GPS recerver 

University of Hawaii 

MeaHuaka 

2006 

N/A 

Center for Robotic Exploration and Space Technolo 

Gene Sat-1 

2006 

N/A 

California Polytechnic Institute 

CP4 

2007 

magnetometer 

Aerospace Corporation 

AeroCube-2 

2007 

N/A 

The Boeing Company 

CSTB-1 

2007 

Sun sertsors. magnetometer 

Tethers Unlimited 

MAST 

2007 

GPS re<»ver 

California Polytechnic Institute 

CP3 

2007 

Magnetometers 

University of Louisiana 

CAPE-1 

2007 

N/A 

University of Sergio Arfaoleda 

Lbertad-1 

2007 

GPS recever 

University of Toronto 

CanX-2 

2008 

Sun sensors, magnetometer 

Tokyo Institute of Techttology 

II 

2008 

Gyro. Magnetometer. Stxi sensor. Earth sensor 

Delft Univefsity of Technology 

Delli-C3 

2008 

Sun sensor 

AHrorg Universitv 

AAUsat-2 

2008 

N/A 

Fachhochschule Aachen 

Compass One 

2008 

Magnetometer, sunsensors 

Nihon University 

SEEDS2 

2008 

Gyro sensor, magnetometer 

Hawk Institute of Space Sciences 

Hawksat-1 

2009 

N/A 

Santa Clara University. NASA Ames 

Pharmasat-I 

2000 

N/A 

California Polytechnic State University 

PolysatCPfl 

2000 

Magnetometers 

Texas A&M University 

Ao9esat-2 

2000 

N/A 

Universitv of Texas at Austin 

BEV01 

2009 

N/A 


Table 3. CubeSat Attitude Determination Methods. After [5]. 
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Another dimension to this is the nearly complete lack of a pre-packaged 
ADS. Until August of 2009, there was only one ADS available on the market. It 
was the Pumpkin IMI-100 ADACS (Attitude Determination and Control System). 
Since then, another system has been introduced. It is the SFL CubeSat 
Compact 3-Axis Attitude Actuator and Sensor Pack with Sinclair Interplanetary 
(SFL ADCS). Both of these will fit the CubeSat Standard and provide attitude 
determination and control of about one degree. They differ in included 
equipment and performance. 

1. Pumpkin IMI ADCS 

Pumpkin, Inc. sells this system as a kit [7]. It includes the IMI-100 ADACS 
developed by IntelliTech Microsystems, Inc. (IMI), a magnetometer, and custom 
software. A separate sun sensor can also be supplied, but the standard setup 
uses only designated solar cells on the CubeSat. This sun sensor setup turns 
the entire CubeSat into a 3-axis sun sensor. The software simply takes the 
calibrated analog signal from the three sides of the CubeSat that can see the sun 
and uses that as the sun vector. The magnetometer is a PNI Micromag3 with a 
32 nT resolution. It must be mounted outside of the IMI-100 to get uncorrupted 
data. The IMI-100 control equipment consists of three reaction wheels and three 
torque coils both orthogonally arranged. The reaction wheels can provide 1.11 
mNms of momentum and 0.635 mNm of torque. A more powerful ADACS is also 
available—the IMI-200, which provides 2.23 mNm of torque and 10.8 mNms of 
momentum [8]. These two systems are compared in Table 4. A microprocessor 
is integrated into the package to perform all the attitude calculations and control 
of the actuators. A pointing accuracy of 1° is advertized with the supplied 
sensors. Higher accuracy is also claimed possible using a Ring Laser Gyro 
(RLG) and a star tracker. 
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Table 4. Specifications of IMI-100 and IMI-200. After [7][8]. 



Figure 4. 2U CubeSat with IMI-100 ADACS (left) and IMI-100 
ADACS alone (right). From [7], 


2. SFL ADCS 

This system is also an actuator and sensor package designed for 
CubeSats [9]. This ADCS has three reaction wheels and three magnetorquers, 
arranged orthogonally, like the IMI. It also comes with a magnetometer and six 
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sun sensors. It, however, does not come with a computer or software. The 
reaction wheels are made by Sinclair Interplanetary [10] and are comparable to 
the IMI-200 with 10 mNms of momentum and 1 mNm of torque. There are four 
sun sensors integrated into the main box with two additional sensors provided to 
be placed on the CubeSat elsewhere. The magnetometer comes on the 
deployable boom to prevent magnetic interference from the satellite. Actuator 
and sensor specifications are listed in Table 5. This ADCS also has an 
advertized pointing accuracy of 1-2° rms. 


SFL 

Reaction Wheel 

Value 

Unit 

Magnetorquer 

Value 

Unit 

ADCS 

Momentum Storage 

10 

mNms 

Dipole 

~0.1 

Am^ 


Max Torque 

1 

mNm 





Dimensions 

5x5x3 

cm 

Dimensions 

8x8x0.4 

cm 


Mass 

120 

g 

Mass 

100 

g 


Power (Max) 

0.7 

W 

Power (Max) 

170 

W 

Sun Sensor 

Value 

Unit 

Magnetometer 

Value 

Unit 


Resolution 

<0.5 

o 

Dynamic range 

±100 

pT 


Accuracy 

<1.5 

o 

Resolution 

15, 30 

nT 


Dimensions 

3x3x1 

cm 

Dimensions 

2x4x2 

cm 


Mass 

<6 

g 

Mass 

<30 

g 


Power 

<0.175 

W 

Power 

<0.025 

W 


Table 5. SFL/Sinclair ADCS Specifications. After [9][10]. 


A note for both of these systems is they only provide good attitude 
determination in sunlight. So, if the CubeSat goes into eclipse, which nearly all 
will, the attitude solution is either degraded, propagated with a model, or simply 
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set to a nominal value. This issue is being addressed now by different 
institutions by incorporating other sensors, i.e., Star Tracker or Horizon Sensor, 
and new algorithms. 

E. TINYSCOPE REQUIREMENTS 

The requirements for TINYSCOPE were originally derived in [11], further 
refined in [12]. These requirements are the basis of this ADS design. However, 
the research here will benefit any nano or pico-satellite that requires either high 
accuracy attitude knowledge or data fusion of several different attitude sensors. 
The primary objective of this thesis is to design and implement a Kalman Filter 
that uses several separately developed techniques to combine different types of 
sensor data in a computationally efficient manor and produce an attitude solution 
that is significantly more accurate than any one of the sensors alone. The 
secondary objectives of this thesis are to implement realistic simulations of the 
various sensors used by the Kalman Filter and to generally increase the fidelity of 
the current Simulink® simulation model. 
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II. SELECTION OF ATTITUDE DETERMINATION HARDWARE 


A. PURPOSE 

The development of high accuracy attitude determination for CubeSat 
class satellites is still new. The accuracy that is demanded by TINYSCOPE is 
greater than any other satellite of comparable size that has flown thus far. This 
means that better sensors must be developed or sensors that existed must be 
used in such a way to obtain better results. 

B. BASIC ADS DESIGN 

The ADS developed for this thesis uses commercially available (COTS) 
sensors. Only minor modification will need to be done to select parts. It will also 
be robust using five different types of sensors to ensure continuous solutions and 
allowances for sensor failure. The key part of this is fusing the data in such a 
way that the attitude knowledge is more accurate than any of the individual 
sensors. The ADS also cannot be dependant on one sensor to provide 
acceptable accuracy. Having a number of different types of sensors and using 
an Extended Kalman filter is the selected solution for the ADS. 

The sensor types selected are mostly inertially referenced. The one 
exception to this is the magnetometer, which is referenced to the Earth’s 
magnetic field. The other sensors will include gyroscopes, sun sensors, and star 
trackers. A GPS receiver will also be incorporated. However, it will not be used 
for the attitude determination; it will be part of the orbit determination system 
(ODS). This system will not be developed in great detail for this thesis. 

C. COMPONENT SELECTION 

An extensive search for sensors that would provide high accuracy attitude 
knowledge for a CubeSat class satellite was conducted over the six-month period 
of September 2008 to March 2009. There were four primary criteria used to 
determine what sensor would be most useful. 
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1. Size. Given the CubeSat standard sensors need to be able to fit inside 
less than 1U. 

2. Power. A CubeSat can only produce approximately twenty watts of 
power total. 

3. Weight. Due to the CubeSat standard. 

4. Performance. This was the deciding factor for many of the sensors. 

Fortunately, the first three criteria tend to go hand in hand, simplifying the 
selection process. Other considerations were cost, the material properties, and 
electrical interfaces of the sensor, but these were lower order factors. 

1. Inertial Measurement Unit 

The Inertial Measurement Unit (IMU) is one of the key components to the 
ADS. It is actually a combination of two different sensors. The IMU has 
gyroscopes that measure the angular rate of the spacecraft and accelerometers 
that measure acceleration. It will provide attitude information from this sensor 
suite in a near continuous fashion. The EKF will specifically use the gyroscope 
information to estimate the attitude when the other sensors are not providing 
information. This will be discussed in greater detail in Chapter VI. The 
accelerometers are not currently being used, but could be integrated into the 
attitude estimation or the orbit estimation later. This type of IMU is used often on 
unmanned aerial vehicles (UAV) which do not generally require highly accurate 
sensors for flight control. Therefore, there is no real driver for the industry to 
produce very accurate, low noise gyroscopes and accelerometers. The 
gyroscopes performance was graded primarily on two characteristics: the 
resolution and the noise. The resolution had to be smaller than the required 
angular rate knowledge. The noise only had to be as small as possible because 
it would be compensated with filtering. 

From the beginning, commercial MEMS sensors were the primary focus of 
the selection process. In general, they are the only gyroscope technology that 
meets the first three criteria for this sensor. The problem is they generally have 
more noise than the other types of gyroscopes do. The IMU that was eventually 
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selected was the Analog Devices ADIS16405 (see Figure 5). A summary of the 
IMU specifications can be seen in Table 6 and the full specification sheet is in the 
Appendix. This IMU has several very useful features that could be exploited in 
the future. It has three dynamic ranges that are set be software. This means 
that during different modes of satellite operation the dynamic range could be 
changed to suite the expected angular rate range. The ADIS16405 also has low- 
pass filtering. A Bartlett window is provided by two cascaded averaging filters. 
The number of taps for each averaging stage can be set through software. Bias 
compensation can be done with either an automatic routine on the IMU or 
through manual settings. Both are done through software settings. 


Gyroscope 

Value 

Unit 

Dynamic Range 

±300, ±150, ±75 

7s 

Scale Factor 

0.05, 0.025, 0.0125 

7s 

Initial Bias Error 

±3 

o 

Bias Stability 

0.007 

7s 

Angular Random Walk 

2 

7Vhr 

3 dB Bandwidth 

330 

Hz 

Dimensions 

23 X 23 X 23 

mm 

Mass 

16 

g 

Power 

0.35 

W 


Table 6. Analog Devices ADIS 16405 Gyroscope Characteristics. 

After [13]. 
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Figure 5. ADIS16405 with Evaluation Board. 

2. Magnetometer 

This sensor provides coarse attitude information to the ADS. 
Magnetometers used for attitude determination alone cannot achieve highly 
accurate results. This is because they need a reference model to compare the 
measurement to produce a usable vector. This magnetic field difference vector 
can then be used to determine the spacecraft’s attitude with respect to the Earth. 
The magnetometers have become very accurate and small; however, the 
reference models are not very accurate. The best models use a twelve pole 
equivalent model, but even these cannot completely describe the Earth’s 
complex and time changing magnetic field. One degree rms is usually 
considered to be about as accurate as a magnetometer ADS can get. It has also 
been shown that information provided by the magnetometer can be used to 
determine the satellite orbit [14], This will not be pursued in this thesis, but could 
be added to the ODS later. 

The magnetometer was graded on noise and resolution like the 
gyroscope. The resolution had to be small enough that the changes in the 
Earth’s magnetic field could be detected in low Earth orbit. The noise had to be 
low enough that it did not drown out the measurements, but it too is filtered. The 
magnetometer selected is actually integrated into the ADIS16405 IMU. It 
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provided comparable performance to other magnetometers considered at no 
extra cost. The characteristics of the magnetometer are summarized in Table 7. 


Magnetometer 

Value 

Unit 

Dynamic Range 

±3.5 

gauss 

Scale factor 

0.5 

mgauss 

Initial Bias Error 

±4 

mgauss 

Output Noise 

1.25 

mgauss rms 

Dimensions 

23 X 23 X 23 

mm 

Mass 

16 

g 

Power 

0.35 

W 


Table 7. Analog Devices ADIS16405 Magnetometer Characteristics. 

After [13]. 

3. Sun Sensor 

The sun sensors provide an important measurement for sun pointing 
mode and coarse measurements during the normal mode of operation. Multiple 
sun sensors will be needed because each one can only provide two axes of 
information. They are typically arranged to provide full (or nearly full) spherical 
coverage. 

Selection for the sun sensor was primarily determined by accuracy and 
size. Most sun sensors that have higher accuracy are much too large for a 
CubeSat. Other sensors that were either very small or integrated into solar cells 
were too inaccurate. The sun sensor selected is the Sinclair Interplanetary SS- 
411 Sun Sensor. The characteristics are summarized in Table 8 and the full 
specification sheet is in the Appendix. 
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Sun Sensor Parameter 

Value 

Unit 

Accuracy 

±0.1 

o 

Field of View 

±70 

o 

Bandwidth 

5 

Hz 

Dimensions 

34 X 32 X 21 

mm 

Mass 

34 

g 

Power 

0.075 

W 


Table 8. Sinclair Interplanetary SS-411 Characteristics. After [15] 



Figure 6. Sinclair Interplanetary SS-411. 
From [15]. 


4. Star Tracker 

The star tracker provides the high accuracy, or fine pointing, information 
for the ADS. These sensors are essentially cameras that take pictures of stars. 
They use the star light intensity and the relative positions to determine what stars 
are in the picture. Once this is established, through one of many possible 
algorithms, the sensor attitude can be determined. Each star tracker generally 
provides three axes of information, so only one unit is necessary. The accuracy 
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of a typical star tracker on a large spacecraft is a few arcseconds (arcsec). This 
accuracy necessarily degrades as the sensor gets smaller and the optics less 
refined. 

The criteria for selection of the star tracker were primarily size and power. 
This is because very few of these sensors can even fit in a CubeSat. Ultimately, 
though, accuracy and availability were the determining factors. A rough estimate 
of the required accuracy of less that 100 arcseci was made during the selection 
process. The selection based on these criteria became easy because at the time 
only one available product existed, the Commtech AeroAstro Miniature Star 
Tracer (MST) [16] (see Figure 7). Unfortunately, even though this is a relatively 
inexpensive star tracker, it is still way outside a university budget, at roughly 
$ 250 , 000 , 002 . Therefore, this sensor will not be purchased, but its 
characteristics are summarized in Table 9, and the full specification sheet is in 
the Appendix. 



Figure 7. Commtech AeroAstro MST. From [16]. 


1 100 arcsec = 0.028° 

2 $250k is a ROM quote, but is (currently) accurate and includes the baffle design. 

19 




Star Tracker 

Value 

Unit 

Accuracy 

±70 (3-axes 3a) 

arcsec 

Sensitivity 

Up to 4*^ magnitude stars 

Dimensions 

5.4x5.4x7.6 

cm 

Mass 

425 

g 

Power 

<2 

W 


Table 9. AeroAstro MST Characteristics. After [16]. 

5. GPS 

The GPS unit will provide position and timing information for the satellite. 
The position information will be the primary part of the ODS. The timing will be 
used for the spacecraft computing capabilities in general. Small size and low 
power were again the main factors in reducing the number of possible 
candidates. Both terrestrial and space qualified units were considered to help 
enlarge the pool of possibilities. The terrestrial units need minor modifications to 
work in space, so this was not of great concern. The GPS receiver selected is 
the NovAtel OEMV-1G-L1-A (seen in Figure 8). It is a terrestrial GPS that has 
been modified and flown in space previously. It was chosen because of its 
relatively low power requirements and good performance characteristics (see 
Table 10). Of interest, the GPS calculations only use approximately 40% of the 
processing power. The version with the Application Programming Interface (API) 
was ordered to allow user defined programs to run on the GPS processor. The 
excess processing power could be used to run the ODS. The full specification 
sheet is in the Appendix. 
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GPS 

1 

Value Unit 

LI SEP3 

1.8 m 

Time Accuracy 

20 ns rms 

Dimensions 

H 

Mass 

46x71x13 mm 

21.5 g 

Power 

J 

1 W 


Table 10. NovAtel OEMV-1G Characteristics. After [17]. 



Figure 8. NovAtel OEMV-1G-L1 with L1\L2 Antenna. 


3 Spherical Error Probability (SEP), means the true position will be within a sphere with the 
given radius fifty percent of the time. 
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III. EXPERIMENTAL CHARACTERIZATION OF SELECTED 

SENSORS 


There are many ways to determine the noise sources of an instrument. 
Two of the most popular are the Power Spectral Density (PSD) and the Allan 
variance. This thesis will use the Allan variance because it is easily computed 
and there is existing data from the manufacturer available. 


A. ALLAN VARIANCE 


The Allan variance, an accepted IEEE standard for gyroscope 
specifications, is a time domain analysis technique that can be used to find the 
characteristics of the noise processes in an instrument. The Allan variance 
technique uses a clustering method. It divides the data into clusters of specific 
length and averages the data in each cluster. It then computes the variance of 
each successive cluster average to form the Allan variance. Each noise source 
has a different correlation time. By choosing the correct correlation time or 
cluster length, the desired noise source variance can be calculated. The Allan 

variance is typically plotted as the root Allan variance, cr^ = on a log-log 


scale. The different noise sources can be discriminated by examining the 
varying slopes of the root Allan variance. A more in depth discussion of Allan 
variance can be found in annex C of [18] or [19]. The general equations for the 
Allan variance forming K clusters from N data points taken at fs samples per 
second with M points per cluster are as follows: 


Cl>, 


yf M 


k = XZ...K 


M 


(^1 1 ’ ■ ■ ■ ’ ’ ®/W+1 ’ ■ ■ ■ ’ ®2/W ’ ■ ■ ■ ’ ^N-M ’ ■ ■ ■ ’ 


k=^ 


k=2 


k^K 


(M) - 0), {M)f^ = (M) - 0), {M)f 


Where <> is the ensemble average and = M / fjs the correlation time. 


( 1 . 1 ) 

( 1 . 2 ) 
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ALLAN VARIANCE 



7" 


Figure 9. Piecewise Representation of Hypothetical Gyro in Allan Variance. 

From [18], 


The following equations have been determined to calculate the variance 
for the Angular Random Walk (ARW) and the Rate Random Walk (RRW) noise 
sources in a gyro and can be found in [18]. The coefficients can be directly 
determined from the root Allan variance in Figure 9. 


'ARW 


( \ 


(1.3) 


'RRW 


(r) = 




3 (1.4) 

where N is the variance coefficient atr = 1along the +1/2 slope and K is the 
variance coefficient at r = 3along the -1/2 slope. These two noise sources are 
the primary noise sources being modeled in Chapter IV. 


The Allan variance algorithm used to estimate noise parameters for the 
gyros was from MATLAB®Central, an open exchange of files for MATLAB® users. 
The specific code, allan v1.71, developed by M. A. Hopcroft [20] uses the basic 
method of Allan variance calculation developed in [19] and has been validated 
using the example data form [21], 
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B. CHARACTERIZING 


1. Gyro Noise 

A test was performed on the I MU to measure the noise characteristics of 
the gyroscopes. The test procedure from section 12.11 of [22] was referenced to 
setup and perform each test. The IMU was sampled with evaluation board using 
a USB interface and software provided by Analog Devices, Inc. The evaluation 
setup was performed per [23]. All tests were conducted at ambient temperature; 
no environmental temperature control equipment was used. To account for this 
the temperature was allowed to stabilize and was also recorded for each test. 
The magnetic field was also not controlled due to lack of equipment to do so. All 
testing was done in a static condition. Several tests were conducted using 
different sample rates, number of samples, and filtering settings (number of 
averaging taps). 

• Gyro Data 1: Data were taken from only the Z axis gyro. It was 

intended to verify the method used here by comparing the calculated graph to 

the graph provided on the specification sheet. The following parameters were 
used: 

• Sample rate: 115 Hz (10 ms delay) 

• Sample points: 2,000,100 (~ 4.8 hrs of data) 

• Gyro Range: 300 7sec 

• Filter taps: 1 (minimum low-pass filtration) 

• Gyro Data 2: Data were taken from only the Z axis gyro. It will 

calculate the noise coefficients for the expected operation mode with no 
filtering. The following parameters were used: 

• Sample rate: 111 Hz (10 ms delay) 

• Sample points 2,000,100 (~4.8 hrs of data) 

• Gyro Range: 75 7sec 

• Filter taps: 16 (minimum low-pass filtration) 

• Gyro Data 3: Data were taken from only the Z axis. It will calculate 

the noise coefficients for the expected operation mode. The following 
parameters were used: 

• Sample rate: 111 Hz (10 ms delay) 
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• Sample points 2,000,100 (~5 hrs of data) 

• Gyro Range: 75 7sec 

• Filter taps: 64 (maximum low-pass filtration) 

2. Magnetometer Noise 

The magnetometer noise was tested in the same way as the gyro. The 
same conditions exist for these tests as the gyro noise tests. The only difference 
is the data was collected from the magnetometers and not the gyros. Several 
tests were conducted using different sample rates and number of samples. 

• Magnetometer Data 1: Data were taken from only the Z axis 

magnetometer. It was intended to find the variance in the measurements. 
The following parameters were used: 

• Sample rate: 115 Hz (10 ms delay) 

• Sample points: 2,000,100 (~4.8 hrs of Data) 

• Magnetometer Data 2: Data were taken from only the Z axis 

magnetometer. It was intended to find the variance in the measurements and 
determine the effect of sample rate. The following parameters were used: 

• Sample rate: 210 Hz (5 ms delay) 

• Sample points: 2,000,100 (-2.65 hrs of Data) 

3. Sun Sensor Noise 

The sun sensor noise is previously characterized by Sinclair 
Interplanetary. Two tests are performed to ensure proper calibration of the 
sensor before it is shipped. The first measures the noise across the entire field 
of view. The test consists of a collimated light source moved across the sensors 
full field of view. Readings are taken to determine the accuracy of the sensor at 
each position of the light source. The dominating noise source is a spatial 
distortion due to variations in the optics and sensor electronics. This noise is 
calibrated down to a 0.1 °rms. The second test measures the time varying noise. 
The collimated light source is held in one position and data is taken to determine 
the noise level. This noise is not significant when compared to the 0.1 °rms 
spatial noise. 


26 



4. 


Star Tracker Noise 


The star tracker was not purchased and so cannot be characterized. For 
all modeling purposes, the star tracker noise is derived from the published 
accuracy. 

C. RESULTS & COMPARISON 

This section will show the results of the above data runs. The calculated 
Allan variance plot and the statistics plots will be compared and explained. 

1. Allan Variance for Gyro 

The first objective for this testing was to verify that the calculations was 
providing accurate data. To do this a data set was taken with parameters to 
closely mimic the same conditions of the manufacturer testing, this was Gyro 
Data 1. The calculated Allan variance was then compared to the Allan variance 
provided on the specification sheet. This comparison can be see in Figure 10 
and Figure 11. 

It is clear that the calculated Allan variance is very similar to the 
manufacturer’s Allan variance. The noise coefficients from the calculated Allan 
variance plot were also checked against the specification sheet. The ARW 
coefficient is 0.028 7Vsec while the specification is 0.033 7Vsec. This is 
reasonably close to say that the difference is due to small variations in each unit. 
The bias stability coefficient was also found for comparison. The coefficient from 
the calculations is 0.0052 7sec while the specification is 0.007 7sec. This again 
seems within reasonably manufacturing variation. It is, therefore, reasonable to 
validate the Allan variance calculations and assume further calculations will 
produce reasonable results. 
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Root Allan Variance f/sec) 


Allan Deviation: Z Gyro 300deg/s 1tap (114.9647 Hz) 
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Figure 10. Calculated Root Allan Variance from Gyro Data 1. 
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Figure 11. ADIS16405 Root Allan Variance. From [13]. 

28 

























































































































































































































The next objective is to calculate the noise sources under parameters that 
more closely resemble the actual operating mode. Gyro Data 2 is a data set 
taken under such conditions with minimal filtering applied. This will help 
determine the effect of the low-pass filtering. Figure 12 shows the root Allan 
variance plot and the noise coefficients. 



Figure 12. Calculated Root Allan Variance form Gyro Data 2. 


The root Allan variance clearly maintains the overall shape seen in the 
previous figures. It can also be seen that, except for ARW, the noise has been 
reduced. The minimal low-pass filtering may be the cause of this. The noise 
reduction does not appear dramatic, but could be significant for improved 
performance of the EKF. 

The last set of data. Gyro Data 3, was taken under operating parameters 
with full low-pass filtering applied. It is expected that the filtering will further 
reduce the noise and improve performance. Smaller values for all the noise 
coefficients will reflect this noise reduction. The root Allan variance plot can be 
seen in Figure 13. 
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Figure 13. Calculated Root Allan Variance from Gyro Data 3. 


The bias stability increased back to the value from Gyro Data 1 and the 
ARW stayed the same from Gyro Data 2. The shape of the graph, however, has 
changed. The RRW shown here is probably inaccurate. The +1/2 slope that is 
normally seen on the Allan variance is either very short or is in the uncertainty of 
the Allan variance calculation. The low-pass filtering appears to have reduced 
the RRW and possibly reduce other noises as well. It also is possible the 
increase in bias stability noise is due to the elimination of RRW. This change in 
noise source levels is not well understood. The -1/2 slope line for the ARW does 
not appear to fit as well as it did on the other Allan variance plots. This could 
mean that the ARW has been reduced slightly, at least for higher frequencies. 

These three root Allan variance plots clearly show the expected outcome. 
The low-pass filtering has significantly reduced the noise of the gyro 
measurements. Table 11 shows a summary of the noise coefficients found from 
these three data sets. 
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ARW (°/Vsec) 

Bias Stability (°/sec) 

RRW (°/Vsec^) 

Gyro Data 1 

0.028 

0.0052 

0.00038 

Gyro Data 2 

0.028 

0.0038 

0.00018 

Gyro Data 3 

0.028 

0.0052 

0.000184 


Table 11. Summary of Noise Coefficients. 


2. Verification of Gyro Model 

The noise coefficients from Gyro Data 2 was used in the Simulink® model 
because of its realistic parameters and high confidence in the results. The 
simulated gyro noise and the tested gyro noise can be seen in Figure 14. The 
mean was subtracted from each data set to better compare the random noise 
levels. The Simulated and real noise are very similar, validating the Simulink® 
model. 

Simulated Gyro Rate Tested Gyro Rate 




Time (sec) Time (sec) 


Figure 14. Simulated Gyro without Bias and Actual Gyro Noise. 


4 This number is questionable due to the poor quality of data. 
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3. Statistical Analysis of Magnetometer 

The formulas and physical relationships have not been developed to use 
the Allan variance for a magnetometer. More mundane methods will be used to 
analyze the magnetometer noise. These will be the mean, to estimate the bias, 
the standard deviation and variance, to estimate the other noises. Two data sets 
at different sample rates were taken from the magnetometer. The statistics were 
calculated and the results are compared. 

Figure 15 shows the statistics from Magnetometer Data 1. The mean (in 
green) is -161.1 mGauss, the 3a (in red) or three times the standard deviation is 
4.57 mGauss, and the variance is 2.32 mGauss. The graph on the left is a plot 
of the Z axis magnetometer over time and the graph on the right is a histogram of 
the same data to confirm the calculated statistics. 


Sampled Data from Z Aws Magnetometer 



X 10® Histogram of Z Axis Magnetometer 


5 



-168 -166 -164 -162 -160 -158 -156 -154 

mGuass 


Figure 15. Statistics Graphs from Magnetometer Data 1. 

Figure 16 shows the statistics from Magnetometer Data 2. The mean (in 
green) is -148.9 mGauss, the 3a (in red) ar three times the standard deviatian is 
4.39 mGauss, and the variance is 2.146 mGauss. The graph an the left is a plat 


32 









































































of the Z axis magnetometer over time and the graph on the right is a histogram of 
the same data to confirm the calculated statistics. 


Sampled Data from Z Axis Magnetometer 



X 10® Histogram of ZAxis Magnetometer 

6p-[-^^-[-[-^^-r 



mGuass 


Figure 16. Statistics Graphs from Magnetometer Data 2. 


The standard deviation and variance were very close for both data sets. 
The means, however, were significantly different. This is believed to be due to 
the test setup of the first data set. The IMU was placed too close to some metal 
washers and the ferrous material corrupted the absolute measurement of the 
magnetic field. This does not invalidate the data because it was constant 
throughout the test and therefore only changed the mean. 

4. Verification of Magnetometer Model 

Figure 17 shows the Simulink® model noise compared to the recorded 
noise. The simulated noise has the same statistics as the two hardware data 
sets. This validates the Simulink® Magnetometer model. 
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Figure 17. 



Simulated Magnetometer Noise and Actual Magnetometer Noise. 
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IV. NUMERICAL MODELING AND SIMULATION 


A. PURPOSE 

This chapter describes the general spacecraft model that was originally 
developed in [11] and has been subsequently further developed for this thesis. 
The notable additions include Earth’s magnetic field model, attitude sensors 
models, and the Extended Kalman Filter described in Chapter VI. The model is 
run from a MATLAB® script file to set initial conditions and plot the results. The 
Simulink® model is described first and focuses on the additions to the model. 
The TINYSCOPE Simulink® Model and MATLAB® script files can all be seen in 
the Appendix. 

B. SIMULINK MODEL 

1. Orbital Propagation 

The standard Euler’s Equation(2.1) were used to create the orbital 
propagation block. The position vector R, the velocity vector V, were outputs for 
use in other blocks in the model. This double integrator system allows for the 
propagation of any orbit by altering the position and velocity initial conditions. 
This flexibility is very useful when analyzing a wide variety of possible orbits. 

R = -4R (2.1) 

r 

Other orbital elements are also computed here. The Beta angle {/ 3 ), 
which is the angle between the sun line and the subsolar point (SSP), and the 
True Anomaly (v), which is the angle between the position vector and periapsis, 
are two of these. They are calculated with Equations (2.2) and (2.3). The non¬ 
orbital elements, latitude and longitude are computed here as well for 
convenience. The eccentricity (e) and inclination (/) of the orbit are predefined 
and constant for the purposes of the simulation. 
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/? = sin ^(sint/sin/cos/-i-cosL/sin/sin^-sinL/cos/sin/cos®) (2.2) 

Where the right ascension of the sun in the ecliptic plane is 
t/=t/o+wf,w =0.985648°/mean solar day, the Earth spin axis tilt with respect to 
the ecliptic plane is/= 23.442°, and the right ascension of the line of nodes 

is(y = (Do + (yf, (y = (-9.9639/?^®cos/)/(l-e^)^°/mean solar day. Where the 

Earth angular radius is/? = sin"' / (R^ + /))). 

v = cos - (2.3) 

y er J 

where, E = -///r)R-(R* V)v)/// and r is the magnitude of the position 

vector of the satellite. 

2. Environmental Effects 

This block models the space environment that the satellite will be 
experiencing. It uses models developed from commonly known relationships. 
There are three aspects of the space environment were modeled; they are: the 
Earth’s Magnetic field, the Sun, and the Earth’s atmospheric density. 

a. Earth’s Magnetic Field Model 

The Earth’s magnetic field can be approximated by a magnetic 
dipole. This model is not very accurate because the true magnetic field has 
many variations due to interactions from within the Earth as well as from outer 
space. The dipole is, however, much easier to calculate and is accurate enough 
for the purposes of these preliminary simulations. Equation(2.4)5 is the vector 
form of the magnetic dipole where m is the Earth’s vector dipole moment, //gis 

the permeability of free space, and r is the satellite position unit vector. Because 
the Earth’s magnetic field is displaced by approximately 11.7° from the Earth’s 

5 The delta function is zero for the magnetic field except at the origin of the dipole. Since the 
spacecraft would never be at the center of the Earth, this term was neglected. 

36 



axis, the dipole moment was rotated into the ECEF frame. The transformation 
matrix used was Equation(2.5). The position vector also had to be rotated from 
ECl to ECEF in Equation(2.6) using the Earth’s rotation 
rate, cp = -7.2921 e-005 rad/sec. 

B = ^^(3(m.r)r-m) + ^m^'(r) (2.4) 

'10 O' 

= 0 cos^ -sin^ (2.5) 

0 sin6' cos6' 
cos^3 sin^3 0 

= -sin^3 cos^3 0 (2.6) 

0 0 1 

Another option tried was the World Magnetic Model 2005 from the 
Aerospace Blockset of MATLAB/Simulink® for its high fidelity. It takes the 
altitude, the latitude, and the longitude of the spacecraft as inputs, all of which 
were derived from the position vector produced by the orbit propagator block. 
The output is the magnetic field vector in Geodetic coordinates. Unfortunately, 
this model is too computationally intensive to be practical in the current 
implementation of the spacecraft model. 

b. Atmospheric Density 

A look-up table using atmospheric density data points from 100 to 
700 km approximates the nominal atmospheric density, . This density is used 

as a starting point for Equation(2.7). This equation takes into account the 
differences in density due to the sun and eclipse. It also accounts for the beta 
angle. The initial density is chosen based on the nominal altitude and the 
expected solar activity. 

(2.7) 
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c. 


Solar Simulation 


The sun is simulated using two well known equations that describe 
the Sun Direction Vector of the body with respect to the orbital frame, 
Equation(2.8), and when the satellite is in eclipse, Equation(2.9). 

cos/? sin V 

Sb=Ceo sin/? (2.8) 

cos/? cos V 


COSV' + 


COSyO 

cos/3 


<0 


(2.9) 


3. Dynamics and Kinematics 

This block uses standard equations to calculate the rigid body dynamics 
as well as the quaternion kinematics. 


a. Dynamics 

The Euler Equation(2.10) is used to calculate the nonlinear 
dynamics of the spacecraft. These equations relate the applied torques to the 
spacecraft angular rates. Only spacecraft rigid body dynamics are accounted for 
and all forces are measured along the three principle axes. 

T = Ucb^+{J,-J,)co,co, ( 2 . 10 ) 


b. Kinematics 

This is the standard quaternion Equation(2.11) that relates the 
angular velocity and the current quaternion to produce the derivative of the 
quaternion. This can then be integrated to produce the propagated quaternion. 

[0 ry, -cOy 

0-2 ° ^2 (2 11 ) 

(73 2 o)y -0)^ 0 0)^ 

_C74j -«y -«z 0 Jl^4_ 
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4. Disturbance Torques 

These blocks calculate the three primary disturbance torques. They make 
heavy use of the environmental blocks outputs discussed earlier. 


a. Gravity Gradient Torque 


This disturbance torque is due to gravity acting on the spacecraft. 
Variations in mass occur throughout any satellite and therefore gravity acts more 
strongly on certain parts of the spacecraft than others. This gradient in strength 
of gravitational force causes a torque on the spacecraft body. The standard 
equations for the gravity gradient disturbance torque are; 


' GG 


3ju 




'Ci' 


0 


where, 

C2 

II 

0 

DO 

0 

0 





1 


( 2 . 12 ) 


b. Aerodynamic Torque 


This disturbance torque is due to the inelastic impacts of molecules 
in the upper atmosphere onto the spacecraft. These impacts transfer momentum 
from the molecule to the spacecraft causing a drag force. The net drag force 
acts on the center of pressure located on the spacecraft face in the direction of 
travel. The equation describing this torque on a spacecraft assuming flat 
surfaces is: 




(2.13) 


where, Cp, = R 


CP/ 


'■CM 


and Vr is a unit vector along relative velocity of 


atmosphere WRT the spacecraft. 


H(v,-h,) = 


0 , 


V . n ,.>0 

V.n,<0 
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c. Solar Torque 

This disturbance torque is very similar to the aerodynamic torque. 
The major difference is that the force acting on the spacecraft is due to the 
impingement of light from the sun. The sun vector calculated in Equation(2.14) is 
used for this calculation. 


n* = -pA ^ i(c„- (s. ■ A, ) A )h. (s. a, ) 


^pi ~~ ^CP/ 


(2.14) 


H(s.n,) = 


where, 

‘‘Cf 

1 , 


0 , 


S, n. >0 
S, n. <0 


5. Sensor Models for Simulations 

This block contains four different types of sensors. A magnetometer, a 
gyroscope, a sun sensor, and a star tracker are simulated with noise sources. 
The goal was to create sensor simulations that would accurately mimic actual 
hardware. 


a. Gyroscope 

The gyroscope takes the angular velocity from the Euler Equations 
and adds noise to create a realistic measured® . Two types of noise are added 
to the signal. Angular Random Walk (ARW) and Bias. The ARW is a zero-mean 
Gaussian random noise with a variance determined from Chapter III. The bias is 
modeled as the integration of white noise with a variance called Rate Random 
Walk, also determined in Chapter III. The integrator is initialized at the initial bias 
of the hardware. These terms can be seen in mathematical model of the gyro 
from [24], Equation(2.15). 

® = ® + + 

P = riu 
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Where(y is the true angular rate, (S is the bias, is gyro scale factor errors, is 
gyro misalignment errors, ; 7 js the ARW, and ; 7 „is the RRW. This equation was 
the basis for the Simulink® model of the gyro seen in Figure 18. 



Integrator 1 

Figure 18. Realistic Simulink® Model for a Gyroscope. 

The ARW, RRW, and initial gyro bias can easily be seen in the 
model. The gyro noise sources are “scaled” by T 3 °®like in [24]. This method 
uses the sample rate,T^, to correct the units of ARW (7Vs) and RRW (7Vs^) to 

7s and 7s^ respectively. The result produces realistic noise in the model. The 
misalignment inaccuracies are added with a gain block in the model. The gain 
is/g^g+G, where the diagonal values of G are the percent error in scale factor 

and the off-diagonal values of G are the percent error of misalignment. The 
gyroscope model also simulates the dynamic range of the hardware with a 
saturation block. These all combined to create a realistic measurement (y)) to 
be used in the simulation. 

b. Magnetometer 

The magnetometer takes the B-field in the body frame, 6^, adds 

realistic noise to simulate the actual magnetometer hardware characteristics. 

Equation(2.16) describes the model used as the basis for the Simulink® model in 
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Figure 19. One noise source is added to6^, the output noise level, of the 
magnetometer. The specification sheet lists the output noise in rms, this is 
equivalent to the la value of the noise. The scale factor,4a> and 
misalignment, error were also added into the model. 

+ ^ (2.16) 



Saturation 

Figure 19. Realistic Simulink® Model for a Magnetometer. 

The magnetometer simulates other physical characteristics of the 
hardware. A saturation block is used to limit the range. The misalignment and 
scale factor inaccuracies are added with a gain block using the same equation as 
the gyro model. These all combined create a realistic measured 6^ to be used in 
the simulation. 


c. Sun Sensor 

The sun sensor does not add a Gaussian random noise to the 
measured signal because, as stated in Chapter III, the dominate noise is special 
in nature. Also note there is no bias for a line-of-sight sensor. The special noise 
was approximated with a modified Bessel function, Equation(2.17). The 
important things to model were the wave like structure of the noise and the rms 
error of 0.1 °. The result of the Bessel function was stored in a look-up table 
indexed by the coordinate of the Sun vector on the simulated sensor detector. 
The small angle approximation allows the 0.1 °rms error to be directly added to 
the measured Sun vector. The other part of this simulation is determining if the 
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Sun is in view of the sensor. First, it is determined if the spacecraft is in eclipse. 
This was described earlier in the Environmental Effects section. Next, the sun 
vector must be changed into the body coordinates, S^. This is done by 


multiplying by Cg^^. Then it is rotated into the individual sensors frame. Two sun 


sensors where simulated in this model, but more could easily be included. The 
two sun sensors frames are set to maximize the time the sun will be seen by at 
least one sensor. Then it can be determined if the Sun is in the sensor’s field-of- 
view. This is done by ensuring the is on the correct side if the spacecraft and 


is in the sensor’s field-of-view as in Equation(2.18). 


I r y ^ 


k=0 


f 


v4y 


k\T(v + k + i) 


cos{0) 


(2.17) 


Where Pis the gamma function, v is a real constant, and r and 6 are the radius 
and angle in polar coordinates respectively. 

S,<-jyyd^, (2.18) 

Whereis the location on the detector and the - denotes a unit vector. 


DCM for 
Sun Sensor 1 



Figure 20. Realistic Simulink® Model for Two Sun Sensors. 
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Figure 21. One Sun Sensor Facing Block. 


d. Star Tracker 

The star tracker model used here is a much simpler one than in 
other simulations. This is because the EKF developed here will assume the star 
tracker is a ’’black box.” That is to say that the EKF will only be a user of the star 
tracker solution and not be an integral part of the star tracker itself^. This will 
allow the EKF to work with any COTS star tracker that might be developed later. 

The star tracker will output a quaternion as many star trackers do, 
however, it will use Euler angles as an input. This allows for easy application of 
the known error of the star tracker to the measured data. A Gaussian random 
noise with a variance of ±70 arcsec is added to each Euler angle. The angles 
are then used to calculate the measured quaternion. 


6 The EKF can be easily modified to be an integral part of a star tracker. The EKF developed 
in [27], which this EKF relies heavily on, does just this. 
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Figure 22. Realistic Simulink® Model for Star Tracker. 

6. Gain Scheduled Quaternion Feedback Controller 

This block was developed in [11] and has not been modified for this thesis. 
No simulations will be done using the controller so it will not be included in this 
description. 

C. MATLAB® CODE 

The Simulink® Model used a number of MATLAB® M-files to both initialize 
the simulation and perform calculations. Files that were not developed for this 
thesis will not be explained here. All the code used for this thesis is in the 
Appendix. 

1. TINYSCOPE Main Script 

The TINYSCOPE Simulink® Model is initialized and run from 
TinyscopeMainScript. This file has all the constants and parameters to run the 
simulation in it. It also calls functions to perform specific calculations. 
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2 . 


Euler to Quaternion 


This script is a function that is called both in TinyscopeMainScript and in 
various blocks of the TINYSCOPE Simulink® Model. It takes the Euler angle 
vector and calculates the 3-2-1 sequence quaternion. 

3. Quat2Euler 

This script is used to convert a quaternion to an Euler angle in the 3-2-1 
sequence. 

4. Calculate 6U Spacecraft 

This script was developed in [11]. It generates the moment of inertia, 
center of gravity, and center of pressure for the TINYSCOPE model. 

5. Plotting Functions 

Several plotting function were developed to aid in presenting the 
simulation results. These functions take in simulation data and generate a series 
of plots that display the data in a usable fashion. Plotting functions developed 
here include: 

• PlotOrbit 

• PlotMeasurements 

• PlotEKFErrors 
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V. KALMAN FILTERING APPROACH TO STATE ESTIMATION 


A. BACKGROUND 

The Kalman Filter is a method to recursively estimate the state vector 
using stochastic processes. The filter finds the optimum solution by minimizing 
the mean square error of the estimated state vector with a system model of the 
plant dynamics and sensor noise. R. E. Kalman first developed this method in 
1960 [25]. Since then it has been used in multiple disciplines ranging from signal 
processing to spacecraft control. It has also been expanded upon and further 
developed many times to now include very good nonlinear state estimators. 

The Kalman Filter uses a two-step process of predicting the state vector 
using a system model and updating the state vector using measurements (see 
Figure 23). This means that only the previous state estimate need be stored until 
the next time step. Thus, the recursive nature of the Kalman filter makes it easily 
implemented on a digital computer. 



Time Update 
(“Predict”) 




Measurement Update 
(“Correct”) 



Figure 23. Kalman Filter “Predict-Correct” Cycle. From [26]. 

This also leads to the advantage that multiple sensors of different types 
can be used to update the estimated state vector. Through an error covariance 
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matrix, the confidence level of each measurement is tracked to appropriately 
weight the update. When a low accuracy measurement is used to update the 
state vector, it can be weighted lower than the predicted model; while a very 
accurate measurement is weighted much more heavily then the predicted model. 
The result is an estimation that when properly implemented can provide more 
accurate state estimation then the direct measurements alone. 


B. DISCRETE EXTENDED KALMAN FILTER 


This section will show the structure of the Kalman Filter. Specifically a 
formulation using discrete nonlinear plant dynamic models with discrete sensors 
will be shown. This formulation is of interest because the implementation of the 
EKF developed will be on a microprocessor. The digital nature of the 
microprocessor lends itself to discrete functions. So, the spacecraft dynamics 
and kinematics are approximated with discrete nonlinear equations and the real 
sensors are discrete digital sensors. The derivations have been done in 
numerous publications but for the purposes of consistency and simplicity all of 
these equations are from [27]. 


To begin, the continuous nonlinear model and measurements are defined 
as: 


x{t) = f(x(t).u{t).t) + G{t)vi{t) 

y(() = h(x((),0 + v(0 


(3.1) 

(3.2) 


Then using the first-order approximation of the nonlinear system dynamics f and 
nonlinear measurement h continuous Riccati equations with Equation(3.3), we 
can approximate the discrete Riccati equation with a Taylor-series expansion for 
exp(Frs) to the second term in Equation(3.4): 


F{x,t) 





0,^/ + FT, 


5h 


dx 


(3.3) 

(3.4) 


Applying Equations (3.3) and (3.4) to Equations (3.1) and (3.2) leads to the 
discrete nonlinear system dynamic and measurement equations: 
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+ (3.5) 

y/<=h,K)+v, (3.6) 

Where Wk & Vk are zero-mean Gaussian noise process with covariances given by 
the expectation equations: 



I 0,k^j 

\o,,k = j 



j 0,k^j 
\R„k=j 



(3.7) 

(3.8) 

(3.9) 


The Qk matrix accounts for state process noise while the Rk matrix accounts for 
expected measurement noise. 


The current state will be propagated by estimating the truth model of 
Equation(3.5) with: 

Xk.i=®kX;+r,u, (3.10) 

The current state will be updated with a measurement, yk, in: 

= X’(3.11) 

The gain Kk changes with time properly weighting the relative confidence of the 
accuracy of the propagated state verses the measured state. To find Kk first the 
state error and error covariance matrixes must be defined: 



(3.12) 

(3.13) 


Substituting Equations (3.5) and (3.10) into Equation(3.12) and substituting the 
resulting equation into Equation(3.13) leads to: 

( 3 . 14 ) 
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Because Wk and are uncorrelated the terms £|w^x^^| = £|x^w[| = 0. To 


find the updated error covariance matrix, substitute Equation(3.6) into 
Equation(3.11). Then substitute the resulting equation into Equation(3.12) and 
reduce, leads to: 






(3.15) 


In order to actually calculate the gain K, the trace of the updated error covariance 
matrix must be minimized. Solving gives: 


Kk-PkHl 


(x»)k(x;)n-H[(x,) + R, 


n-1 


(3.16) 


The Continuous-Discrete Extended Kalman Filter is summarized in Table 
12. The initialization of the EKF is important because the estimated state is 
assumed to be close to the true state. Large initial condition errors on some 
nonlinear plants can cause instabilities in the EKF. This can be avoided through 
simulation and testing. 


Discrete Extended Kalman Filter 

Model 

=®A+r,u,+Y,w, 
Yk=Hk{^k) + ^k 

Initialize 

x('o) = Xo 

P„=£jx(1„)x^((„)) 

Gain 

K,=P;Hl(K) 

H,(x-)PX(x;) + R.]'’ 

Update 

x; =x 

p^ - 

/;+K[y»-h(x;)_ 

Propagation 

p,;, =+ tfijl 


Table 12. Discrete Extended Kalman Filter. From [27]. 
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C. CHALLENGES OF MULTIPLE SENSOR SYSTEM 

In the formulation of the above EKF, the sensors are modeled as discrete 
measurements. This is a very good model for digital sensors or sampled analog 
sensors. This type of sensor is ideal to be used with a Kalman filter implemented 
on the microprocessor. There is, however, an inherent problem. Different types 
of sensors produce solutions at different rates. This causes a problem for a 
Kalman filter that has to have all of the measurements available to update the 
state vector. A solution to this is to use superposition for the updates. This is 
possible for an EKF because although it is modeling a nonlinear system it 
linearizes the propagation and update equations about the current state estimate. 
This technique of using superposition, was first suggested by James Murrell in 
[28]. It has since been applied many times and in different ways. 

The essence of the technique is to update the gain, error covariance, and 
state error vector with each successively available measurement. Once all of the 
measurements have been taken into account, the EKF will propagate the 
estimated state and covariance matrix until the next measurement or set of 
measurements are available. Interestingly, this also greatly reduces the 
computational burden. Instead of calculating a gain matrix that requires an 
inverse of a 3n x 3n matrix, only a 3 x 3 matrix inverse is required n times with 
Murrell’s version. 
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VI. IMPLEMENTATION OF EXTENDED KALMAN FILTER FOR 

MULTI-RATE SENSORS 


The Discrete Extended Kalman Filter discussed in Chapter IV will form the 
bases of the EKF developed for this thesis. The EKF developed here will use 
ideas and equations from several sources. These include [14], [27],[29], and 
[30]. The combination of these different sources helps create an EKF that can 
handle data from gyros, star trackers, sun sensors, and/or magnetometers in an 
efficient way. The goal of this EKF is to produce an accurate attitude estimates 
using a gyro and any combination of other sensors. It will also be 
computationally efficient for implementation on a low power microprocessor. 

A. MUTIPLICATIVE QUATERNION EXTENDED KALMAN FILTER 

Because this EKF will use different types of sensors running at different 
rates and it will be implemented with on-board computation, a version of Murrell’s 
form briefly described in Chapter IV was used. The basic structure and much of 
the derivation of this EKF comes from [27]. The significant changes occur in the 
formulation of the observation matrix Hk and the calculation of the residual, ^, for 
each of the sensor types. To begin the derivation the state vector, Equation(4.1), 
was chosen to be composed of four quaternion elements and three gyro bias 
elements. 

X = [Qi Qs Qs Qa P. Py (4-1) 

Where the quaternion is defined by q = [esin(a/2) cos(a/2)]^ = [g q 4 ]^and 

follows the normalization q^q = 1. The normalization constraint prevents a simple 
calculation of the quaternion error by subtraction. A different approach called the 
multiplicative error quaternion must be calculated. This is defined as 

<^q = q0q' (4.2) 

Where 5c\ = [5g and the inverse quaternion is q"' = [-g q^f . Taking the 


time derivative yields: 
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^q = q(S)q Vq®q ^ 

This eventually gives the estimated quaternion kinematics'^ 


^ = |s(q)« = |Q(«)q 


-(q)- 


^4^3x3 [$" 




Q((y) = 


-[^yx] 


-CO 


CO 

0 


(4.3) 

(4.4) 

(4.5) 


Now to find the discrete propagation of the quaternion kinematics 
Equation(4.4) will be expanded using the power series approach in 
Equation(4.6). 


n{a)t 





{2k)l 



2k+^ 

(2/f + lj 

! 


(4.6) 


Next substituting the identities Q^''((y) = (-1)^||(y||^^/4^4 and 
q 2 /<+i(^)^(- 1 )'<||^|| 2 ''q(^) into Equation(4.6) gives 


n(a)t 


= ' 4 rtZ- 


(- 1 )‘ 


2II II 


2k 


■+ ky 


00 


(- 1 / 


2II II 


-l2k+1 


k^c . . ^ (2/C + 1)! 

Then simplifying with the cosine and sine functions leads to Equation(4.8). 

^1 ^ 


(4.7) 


n(®)f 


2 _ 


UxA COS 






Sin 


+ 0(^(5^- 


mt 


j 


ky 


Now the quaternion propagation can be simply written as 

q/<.i=^(4")q/: 


(4.8) 


(4.9) 


Where ®^and q^^are the post-update estimates and is defined in 

Equations (4.10) and (4.11). 


^ See [27] for more details. 
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n(4*)- 


^ 3 x 3 cos - At -[‘F 


cos -\\d)t\\At 




sin d)t At a. 


(4.10) 


(4.11) 


The discrete sample time is denoted by At. The next propagation models to be 
defined are the angular velocity and the gyro bias. These models follow the 
standard EKF formulation given the post-update bias . 

= A=A (4-12) 

The error covariance propagation follows the same approach as in Chapter V. 
The discrete model is derived using a power series like was done for the 
quaternion update. This leads to Equation(4.13) whereandare 
defined in Equations (4.14), (4.15), and (4.16). 


^k+^ ~ + ^k^k^k 

r®ii ^ 12 ] 


(4.13) 


^k = 


®21 ®22 


®12 = [^k X_ 


sin( 

\ nA x1 ^ 
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yojk XJ 

K1 
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0 

(/) 

^ + 

> 

*^-4x3^ 

2 

l^/c 


0)21=1 


0)22 = 


■[(Uk xj 


1-COs(|«;|^ |Afj 


d)^ Af-sin^|(5;^| Afj 


(4.14) 


4x3 ^3x3 

^3x3 4x3 


(jAt + —(7^At 4^3 -\—(7,At\l. 


- -^<Af 4,3 

V / 


'3x3 


2 


(4.15) 


(4.16) 
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Where and cr^ are the variances 

Now the update equations need to be determined. To start, the estimated 
error state definition isAx^s . Using the small angle 

approximation 2 and q 4 ~^, the four state quaternion has been replaced 

by the three state Euler error angle vector. This minimizes use of the factors 
and 2 in the EKF. It also allows the 3a bounds to be directly calculated from the 
error covariance matrix. 

The equations for the calculation of the observation matrix, Hk, for the 
different sensor types will now be derived. First, it must be determined what 
each sensor is actually measuring. The star tracker will provide a quaternion 
solution as its measurements, each sun sensor will provide a body vector 
referenced to the sun, and the magnetometer will also provide a body vector but 
referenced to the Earth. If the star tracker were to provide body vectors 
referenced to the stars, its update would be the same as the sun sensor update. 
This dose, however, require access to the star catalog and since the star tracker 
is being treated as a black box this access cannot be assumed. Next, the 
measurement model will be developed. The definition of the measurement is 

y/< + (4.17) 

For the star tracker that provides a quaternion as the measurement, 
Equation(4.17) becomes 

(4.18) 

The actual attitude quaternion, qk, is related to the propagated measurement, q^ , 
through 

q, =Aq,0q, =[s(q,) q,]Aq, (4.19) 

Using Equation(4.19) in the observation matrix definition Equation(4.20) leads to 
the observation matrix for the star tracker becoming 

(4.20) 
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(4.21) 


^S(q,) 03,3 

The sun sensor provides body vectors for which the true and estimated vectors 
are defined as 

b../((q.)r b;=/((q;)r (4.22) 

Where the actual attitude matrix, >4 (q^,), and the estimated attitude matrix, , 

are related through 

/t(q.) = /t(<yq,)/t(4-) (4.23) 

With the error-attitude matrix approximated with 

^(<5q«)“'3x3-[a^ (4.24) 

By substituting Equations (4.23) and (4.24) into Equation(4.22) and using the 
result in Equation(4.20) the observation matrix for the sun sensors becomes 

"H(x;) = [[4(q.-)rx] 03 , 3 ] (4.25) 

Finally, the magnetometer also measures a body vector as defined in 
Equation(4.22). However, due of the highly nonlinear nature of the Earth’s 
magnetic field and the resulting uncertainty in the estimated magnetic field it will 
be developed using the perturbation technique. Start with the equations for 
propagated and measured magnetic field respectively. 

e(x() = D;e,+AD;e,+v, e„ = D;e,+v„ (4.26) 

Where D'is the transformation from inertial to body coordinates, AD'S, is the 
estimation error, and and are the errors in the magnetic field model and 
measurement respectively. Defining v = - Vp, the residual becomes 

z = -AD'B,+v (4.27) 

The transformation error, AD', can be defined as the difference between the true 
inertial to body transform and the estimated inertial to body transform. The 
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estimated transform can be written D^^D'. Where D^\s the transform from the 
body to computed coordinates. Now, assuming Df’ is made of small angles, the 
transformation error can be rewritten as 

ad;=(/-|;a«; x])d; -d; = -[a«; x]d; (4.28) 

Where [Aa^x^is the cross product matrix of Euler error angles. Substituting 
Equation(4.28) into Equation(4.27) leads to 

z = [[-B,x] 03 , 3 ]Ax;+v (4.29) 

The observation matrix follows in Equation(4.30), since 6 b is unknown 6 m is 
substituted. 


H=[[-e„x] 03 , 3 ] (4.30) 

Finally, the remaining update equations need to be determined. The 
standard EKF form for the error state vector update is . 


= K, 




(4.31) 


Where, yk is the measurement and h^^x^^is the estimated measurement. For 
the star tracker the estimated measurement is the propagated quaternion, . 
For a body vector it has a form (x^) = , where r is the propagated 

vector to the reference body, i.e. the Sun. The bias update is 

(4.32) 

The quaternion update must us quaternion multiplication like the quaternion error 
did in Equation(4.2). The quaternion update equation is 

(4.33) 


The error covariance matrix update and the gain calculation, both have the same 
form from Chapter V. 


= 

'k 




K,=P,-Hl (xi)[H,(x;)P,Hl(x,) + R, 
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(4.34) 

(4.35) 



The Discrete Multiplicative Quaternion Extended Kalman Filter equations 
are summarized in Table 13. 


Discrete Multiplicative Quaternion Extended Kalman Filter 

Initialize 

q(*o) = qo. P(K) = Pa 

P(K)-P« 

Propagate 

h = k 

^k+^ ~ + ^k^k^k 

Gain 

K,=P,Hl(i;)[h 

H=[[ 

*^3x3 

iA(q,)rx] 03,3 

~^m 03^1.3 J 

Update 

Ax>_ 

Ax; = K, 

k- 

yk-K(^^k)_ 

A + 

% 9/< 

- [/- 

'k ' 

“T OCZi. 

2 " 


Table 13. Discrete Multiplicative Quaternion Extended Kalman Filter. 

After [27]. 
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B. IMPLEMENTATION 


The EKF described in Table 13 is incomplete without some of the 
particular implementation items that must be addressed. The initialization, 
expected measurement noise, and quaternion normalization are a few. all be 
covered to help complete the full development of this EKF. 

1 . Initialization 

In an EKF initialization can be very crucial to the performance. Flere the 
initial state was chosen to beXo=[qQ where 1]. The error 

covariance matrix is a little more involved. Since all of the noise sources in the 
sensors and propagation model have estimated variances (from Chapter III), an 
estimation can be calculated for the error covariance. Using Farrenkopfs 
steady-state analysis from section 7.2.4 in [27], one can calculate a rough 
estimate to the error covariance. The following equation was used to estimate 
Po. 

Where cr„ = sjo-lj + als + o-h^g , cr^is the variance associated with;/^,, and cr^ is the 
variance associated with //^both from Equation(2.15). 

2 . Measurement Noise 

The expected measurement noise was calculated for each sensor in 
Chapter III. They are used in the R matrix when the Kalman gain is calculated. 
Each estimated measurement noise matrix follows the form 

(4.37) 

The actual values of noise variance used are summarized in Table 14. 
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Sensor 

Value 

Units 

Star tracker 

O-,", =1.2797E-8 

rad 

Sun Sensor 

^,"3 =3.0462£-6 

rad 

Magnetometer 

al^=^.25E-6 

tesla 


Table 14. Measurement Noise Variance Values. 


3 . Quaternion Normalization 


The quaternion update used in the formulation of the EKF, Equation(4.33), 
is only guaranteed to be a unit vector to within the first-order. This means 
normalization needs to occur to prevent the build up of computational errors. 
Two methods were used to help maintain a normalized quaternion. The first 
method used does not formally normalize the quaternion, but reduces the error to 
order/32. This factor is suggested in [29], computed in Equation(4.38), and is 
multiplied by the quaternion update when falls between two predetermined 


error bounds. 


1 + 3q,^q 


k y 


(4.38) 


The second method is the brut force normalization per Equation(4.39). 
The full normalization is only done when q^^^q^, becomes greater than the upper 


error bound. This reduces the computational burden for on-board processing 
because it is simple a scalar multiplied by a vector the majority of the time. How 
often this is done should be determined through simulation to establish an 
acceptable error tolerance for the application. 





(4.39) 
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4. 


Murrell’s Version 


The application of superposition to the EKF was discussed in Chapter IV. 
To actually implement it in the EKF some structural changes were made from the 
standard form. The most significant change, however, is in the estimated error 
state update in Equation(4.31). The update now must include any previous 
estimates of the error state from that set of measurements and only add error 
unaccounted for in these previous estimates. Equation(4.40) replaces 
Equation(4.31) in the above EKF. 

Ail = + K, [y, - h, (x,) - H,Ai ,] (4.40) 
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VII. SIMULATION RESULTS 


The simulations performed here are only preliminary. They have not been 
verified through a Monte Carlo simulation or hardware-in-the-loop simulations. 
No disturbance torques are applied during any simulation. No control torque is 
applied during any of the simulations, as only the error between the true value 
and the estimated value is needed. This means the no particular attitude is 
actively maintained, however, the satellite will be initialized nadir pointing and will 
rotate at orbital velocity. This should maintain a nadir pointing attitude because 
there are no disturbances. During each simulation selected sensors will be used 
in the EKF to estimate the attitude. The estimate is compared to the actual 
attitude in Figure 24 and Figure 25 to form the errors. Table 15 shows the 
simulation conditions for each run. Additional Figures of simulation results are in 
the Appendix. 


Simulation 

Sensors* 

(ST, SSI, SS2, 

Mag) 

Initial 

(inertial) 

Initial 

Attitude 

(inertial) 

EKFio 

(Hz) 

Sensorii 

(Hz) 

1 

(1.0,0,0) 

( 0 , 27 t/P, 0 ) 

(0.0,0) 

20 

1 

2 

(1.1.1.0) 

( 0 , 27 t/P, 0 ) 

(0,0,0) 

20 

1 . 5,5 

3 

(0.1.1.1) 

( 0 , 27 t/P, 0 ) 

(0,0,0) 

20 

5 , 5,20 

4 

(1.1.1.1) 

( 0 , 27 t/P, 0 ) 

(0,0,0) 

20 

1 , 5 , 5,20 


Table 15. Simulation Conditions. 


8l=On,0=Off 

9 2tt/P is the orbital velocity wrt the inertial frame. 

The simulation and the EKF ran at the same speed. 

11 The sample rates of each selected sensor is shown in the order of column 1. 
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Figure 24. True Euler Angles in Inertial Framei2. 


X True Gyro Rate 



Figure 25. True Angular Rates in Inertial Frame. 


12 The jumps in angle is due to the range limitations of the sine and tan2 functions. 
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A. EKF PERFORMANCE WITH NOISY STAR TRACKER 

During this simulation, the EKF will only use the quaternion measurements 
from the star tracker. This will set a baseline for other configurations to compare 
against. The expected Euler angle error is less than 0.026 ° 1a. The rms value 
of the error will be calculated and used for comparison. The 3a baundaries (red) 
are directly calculated frcm the errcr covariance matrix. Figure 27 and Figure 28 
show that the Euler error and the bias error are in fact well bound. They also 
show that the rms error is within the expected limits. An unexpected result is 
also apparent. A spike in the error covariance, along with degradation in the 
Euler error at -2853 sec occurs. It has been determined the cause of this error is 
due to the scalar quaternion equaling zero at this point (see Figure 30). This 
causes the gain to become very large and in turn cause the increase in P. Figure 
29 and Figure 31 show a weakness of the EKF, it does not filter the gyro noise 
well because it is used in the propagation process and not as a measurement. 
This accounts for the minimal improvement in the Euler angle estimate over the 
star tracker expected error. 
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Figure 26. Simulation 1 Euler, Bias, and Angular Rate Error. 
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Figure 27. Simulation 1 Euler Angle Error with 3a Boundaries. 
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Figure 28. Simulation 1 Gyro Bias Error with 3a Boundaries. 
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Figure 29. Simulation 1 Gyro Rate Error. 

Figure 30 shows how the Gaussian noise in the Euler angles is no longer 
Gaussian in the quaternion. Figure 31 can be used for comparison with the rate 
error in Figure 29. 
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Simulation 1 Star Tracker Quaternion Measurements. 
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Figure 30. 
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Figure 31. Simulation 1 Gyro Rate Measurements. 

B. EKF PERFORMANCE WITH STAR TRACKER AND SUN SENSORS 

This second simulation makes the sun sensor measurements available to 
the EKF as well as the star tracker quaternion. It is expected that an additional 
measurement will improve the estimate. Again the Euler angle error and the bias 
error are well bound by the 3a boundaries in Figure 33 and Figure 34. The gyro 
bias rms error has improved from simulation 1. The lack of improvement is due 
to the anomaly discussed earlier. The sun sensors are unable to correct for this 
error because the satellite is in eclipse at the time. This can be seen in Figure 36 
and Figure 37, i.e., at -1800 seconds all measurements go to zero. 
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Figure 32. Simulation 2 Euler, Bias, and Angular Rate Error. 
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Figure 33. Simulation 2 Euler Error with 3a Boundaries. 
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Figure 34. Simulation 2 Gyro Bias Error with 3a Boundaries. 
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Figure 35. Simulation 2 Gyro Rate Error. 
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Figure 36. Simulation 2 Sun Sensor #1 Unit Vector Measurement. 
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Figure 37. Simulation 2 Sun Sensor #2 Unit Vector Measurement. 
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c. 


EKF PERFORMANCE WITH SUN SENSORS AND MAGNETOMETER 


During this simulation the sun sensor measurements and the 
magnetometer measurements will be used by the EKF to produce the estimate. 
The star tracker quaternion will not be used. It is expected that this configuration 
will take longer to reach steady state and be less accurate than the star tracker. 
Figure 39 and Figure 40 show that the rms errors are significantly less accurate 
than Simulation 1 or 2. What is interesting is that this simulation was more 
accurate than both ADCS “in a box” discussed in Chapter I. It is clear where the 
eclipse occurs on all the Figures due to the increase in errors, however, one 
should note that an rms error of approximately 1 ° rms is still maintained with only 
the magnetometer measurements. 
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Figure 38. Simulation 3 Euler, Bias, and Angular Rate Error. 
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Figure 39. Simulation 3 Euler Error with 3a Boundaries. 
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Figure 40. Simulation 3 Gyro Bias Error with 3o Boundaries. 
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Figure 41. Simulation 3 Magnetometer Measurement. 


Note on Figure 41 that the range has been changed to ±5 mGauss (5x10''^ 
Tesla) for the simulated magnetometer. The hardware range of ±3.5 mGauss 
proved too small and saturation occurred. This caused large errors in the state 
estimation. It was determined that the simulation should be run simulating a 
similar magnetometer with a larger dynamic range. 
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Figure 42. Simulation 3 Error Covariance Matrix Normal. 
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D. EFK PERFORMANCE WITH ALL SENSORS 

Now all of the available sensor measurements are used in the EKF. It is 
expected that this configuration will produce the best estimate. It is obvious, 
however, that this is not the case. The Euler angle error is twice that of 
simulation 1. Figure 44 shows that the Euler error is not very well bound by the 
calculated 3a boundaries. Figure 46 shows that the star tracker quaternion 
anomaly is still present, but no longer dominant. It appears that the 
magnetometer error is dominating the other sensors. The eclipse can also be 
clearly seen from 1800 to 3900 sec. The relatively small error covariance values 
suggest a noise value is not being accurately estimated. This is most likely from 
the magnetometer. It has the highest sample rate and the most error in attitude 
estimation. More work needs to be done on both the Earth’s magnetic field 
model and the magnetometer modeling to improve the performance of the EKF. 
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Figure 43. Simulation 4 Euler, Bias, and Angular Rate Error. 
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Figure 44. Simulation 4 Euler Error with 3a Boundaries. 
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Figure 45. Simulation 3 Gyro Bias Error with 3a Boundaries. 
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Figure 46. Simulation 4 Error Covariance Matrix Normal. 
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VIII. CONCLUSION 


A. SUMMARY 

The results from the simulations clearly show that the EKF developed 
requires more tuning to achieve the expected performance. It has demonstrated 
utility in fusing the various sensors, but has not yet shown a great improvement 
in performance over the star tracker alone. 

The anomaly, due to the star tracker quaternion causing the error 
covariance to dramatically increase, was an interesting result. A solution to this 
problem would be to use the centroids and star catalog of the star tracker like is 
done in [27]. This would avoid the problem encountered when the scalar 
quaternion equals zero. This sets of a chain of events which results in the error 
covariance increasing as well as actual reduced accuracy of the estimate. This 
only seems to occur because the scalar quaternion makes up the diagonal 

component ofS^q^j). That makes the diagonal component of the Hk matrix zero, 

which is then used in the calculation of the gain. The gain then affects both the 
estimated state update and the error covariance update. 

The solution to this problem requires the ADS have full access to all of the 
star tracker calculations and databases. This would essentially integrate the star 
tracker into the ADCS computer. An integrated ADS could then consist of star 
trackers, sun sensors, and magnetometers, as well as other sensors like horizon 
sensors, all connected and run by one processor to produce an attitude estimate. 

The gyro model has been a great success, however, for the ADIS16405 
hardware, modeling the bias instability may be more appropriate then the RRW. 
Modeling the bias instability would accurately account for the variations in the 
bias under all conditions. As was seen in Chapter II, the RRW is not a 
dominating noise in this gyro and may even be able to be filtered out. The bias 
instability, however, remains even with maximum filtering applied. 
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The magnetometer included in the ADIS16405 proved to have too small of 
a range for this application. The placement of the magnetometer is also 
problematic. The IMU is normally situated inside the spacecraft and given the 
small size of a CubeSat; it would consequently be close to magnetic interference 
from batteries and other electronics. Therefore, it is recommended to use a 
separate, external magnetometer in the final design of this ADS. 

The general Simulink® simulation model has been significantly improved 
with the generalization done in the orbit propagation. The addition of the Earth’s 
magnetic field and the improvement of the atmospheric density both provide a 
higher fidelity model of the space environment. 

This thesis has not achieved all that was set out to do, but it has 
accomplished much. It has proven that data fusion is possible using Murrel’s 
technique. It has also set ground work for the continued improvement and 
development of the TINYSCOPE Simulation. Most importantly, it has shown that 
a CubeSat can indeed have high accuracy attitude determination using currently 
existing technology. 

B. FUTURE WORK 

1. Verification and Testing 

This thesis only presents preliminary results. A full Monte Carlo 
verification of the EKF performance should be performed to provide detailed 
performance parameters for the EKF and the simulation. Further verification of 
the EKF design can be performed with a 3-Axis Simulator. Running the EKF in a 
hardware-in-the-loop simulation will provide more realistic performance 
information on the EKF. 

2. Further Develop Simulation 

The sensor models developed for this thesis went a long way to providing 
realistic measurements for the ADCS. Further refinements to the Simulink® 
model can be made however. The star tracker model is still a fairly simple model 
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and can be improved in many ways. Different sensors should also be developed 
such as a horizon sensor. The Earth’s magnetic field certainly needs to have a 
more accurate model developed. A high order IGRF model should be developed 
that will provide a more realistic approximation of the actual magnetic field and 
can still run efficiently in the simulation. 

3. Hardware 

The development of new low cost, low power, miniature attitude sensors is 
accelerating. Continued evaluation of these new sensors should be performed 
as well as testing of these new sensors for possible performance improvements. 
This would include the currently available prepackaged systems for CubeSats 
discussed in Chapter I. The ultimate goal could be to develop a prepackaged 
ADS at NFS that could work with the proposed ACS developed in [31]. This 
would necessarily require implementation of the EKF developed in this thesis on 
a microprocessor like the MSP430 from Tl. 
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APPENDIX 


A. ADDITIONAL SIMULATION RESULTS 


1. Simulation 1 
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Figure 47. Simulation 1 Gyro Bias. 
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Figure 48. Simulation 1 Error Covariance Matrix Normal. 
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Simulation 2 
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Figure 49. Simulation 2 Simulated Star Tracker Quaternion. 
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Figure 50. Simulation 2 Error Covariance Matrix Normal. 
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Figure 51. Simulation 3 Gyro Rate Error. 
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Figure 52. Simulation 3 Sun Sensor #1 Measurement. 
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Figure 53. Simulation 3 Sun Sensor #2 Measurement. 
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Figure 54. Simulation 4 Gyro Rate Error. 
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Figure 55. Simulation 4 Star tracker Quaternion Measurements. 
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Figure 56. Simulation 4 Sun Sensor #1 Measurement. 
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Figure 57. Simulation 4 Sun Sensor #2 Measurement. 
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Figure 58. Simulation 4 Magnetometer Measurement. 
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B. SENSOR DATA SHEETS 


This appendix contains all the data sheets for the selected sensors. 


1. Sinclair Interplanetary SS-411 


SS-411 Two-Axis Digital Sun Sensors 



Accuracy 

±0.1° over+70° Field-of-View 

Earth Albedo Error 

Rejected by internal digital processor 

Bandwidth 

5 vector solutions / second 

Command / Telemetry 

UART. SPI, FC, CAN or RS-485 

Mechanical 

34 mm x 32 mm x 21 mm, 34 g mass 

Outer Surfaces 

Low emissivity gold, scratch-proof sapphire 

Supply Voltage 

5.0 to 50.0 V 

Supply Current 

7.5 mA avg, 26 mA peak 

Environment 

-25°C to +50°C operating 

16 grms random vibration, 5000 g shock 

20 krad total dose at components 

Heritage 

14 Sensors on-orbit on 4 LEO satellites 

Flight units delivered or built for 20 more satellites 


Sinclair 

Interplanetary 


Contact 
(647) 286-3761 (voice) 
(775) 860-5428 (fax) 
dns@sinclairinterplanetary com 
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AeroAstro Mini Star Tracker 



Miniature Star Tracker 



MST - Horizontal View 



MST - Vertical View 


AeroAstro’s Miniature Star Tracker (MST) design 
is a small, low power (< 2 Watt) star tracker, with 
an accuracy of better than ±70 arc-seconds in all 
three axes (3a). It achieves reasonable star 
tracking accuracy with low mass and power 
consumption at less than half the cost of other 
star trackers. 


The MST is available in a low-cost commercial-off- 
the-shelf (COTS) version tolerant up to 30 krad and 
can be modified for higher radiation tolerance. 
The MST also provides a lost-in-space capability 
and is currently being enhanced to achieve fast 
angular rate sensing. 

The star tracker features a user-definable star 
catalog and powerful hybrid processor. With a 1 
Mpixel CMOS array, the star tracker is sensitive 
up to 4”' magnitude stars. 

Images can be downloaded for ground processing 
and custom code can be incorporated. Built-in 
test includes the ability to upload images and 
verify star tracker performance. 


Making Space For Everyone'" 
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Specifications _ 

Accuracy: Better than ±70 arc seconds, 3-axes (3o) 

CMOS Imager: - 1024 x 1280 pixel array each pixel - 7 
pm square 

Sensitivity: Up to 4th magnitude stars 

Maximum Pitch/Yaw Rate: lOVsec (goal) 

Update Rate: - 1 Hz 

Output: Quaternion, Centroids, and custom 

Stars Tracked: Up to 9 simultaneously 

Star Catalog: 600 and can employ user-defined 
catalogs 

Image Rate: 0 to 24 fps 

Power: < 2 Watts 

Radiation Tolerance: Up to 30 krad(Si), more with 
shielding 

Dimensions: 2”x 2”x 3”(5.4 cm x 5.4 cm x 7.6 cm) 

Mass: 425 g (not including baffle) 

Self Test: Images can be up and down loaded for 
verification 


<:''<} AeroAstro 

20145 Ashbrook Place 
Ashbum, VA 20147 USA 
info^eroastro.com 
vwvw.aeroa5tro.com 

t: 703.723.9800 
f: 703.723.9850 

ftev 07*02 
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NovAtel OEMV-1G 


OEMV-I & 
OEMV‘-1G 


NovAter's OEMV-1 products are high performance, single frequency GPS engines 
featuring additional signals and API support. The OEMV-1 offers integrated LI GPS and 
OmniSTAR and CDGPS while the OEMV-IG offers LI GPS+GLONASS, both in a compact 
size with low power consumption. 


Features 
Ll GPS+Integrated 
OmniSTAR and CDGPS 
on OEMV-1 




The Application Programming Interface (API) functionality introduced on 
NovAtet’s OEM4 receivers is also available on all OEMV receivers. With 
the API library and a standard C/C-f development environment, an 
application can be developed to run directly from the receiver 
platform, eliminating system hardware, reducing development 
Noitel time and resulting in faster time to market. 
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OEMV^-I/IG 


Performance' 


Physical & Electrical 


Channel Configuration 

OEMV-1 OEMV-IG 

14GPSL1 14GPSL1 

IL-band 12GLOL1 

2SBAS 2SBAS 

Horizontal Position Accuracy (RMS) 



OEMV-1 

OEMV-IG 

LI 

1.8 m 

1.8 m 

SBAS' 

0.6 m 

0.6m 

CDGPS' 

0.6 m 

na 

OmniSTARVBS- 

0.7 m 

na 

DGPS 

0.45 m 

0.45 m 

RT-20^“ * 

0.2 m 

a2m 

Measurement Precision 


LI C/A Code 


4 cm RMS 

LI Carrier Phase 

0.50 mm RMS 

Data Rate^ 

Measurements 


50 Hz 

Position 


50 Hz 

Time to First Fix 
Cold Start^ 


60s 

Hot Start* 


35 s 

Signal Reacquisition 

LI 

0.5 s (typical) 

Time Accuracy ’ 


20 ns RMS 

Velocity Accuracy 

0.03 m/s RMS 

Dynamics 

Velocity® 


515 m/s 


Size 46x71x13 mm 

Weight 21.5 g 

Power 

Input Voltage + 3.3 +5%/-3% VDC 
Power Consumption 1.0W (typical) 

Antenna LNA Power Output 
Output Voltage 5 V nominal 

Maximum Current 100 mA 


Communication Ports 

• 2 LV'TTLserial port capableof3(X>to 
921.600 bps 

• 1 L.V-TTLserialportcapableof300to 
230/400 bps 

• 2CANBus^serialportcapableof1 Mbps 

• 1 USB port capable of 5 Mbps 


Input/Output Connectors 

Main 20'pindual row male header 

Antenna Input MCX female 

Environmental 

Temperature 

Operating -40*Cto ^BS^C 

Storage •40*Cto+85*C 

Humidity 95% non-condensing 


Regulatory 

RandomVibe RTCADO-160D(4g) 

Bump/Shock MIL-STD810F (40g) 


Optional Accessories 

• GPS-700 series antennas 

• AIVT-500 series antennas 

• RF Cables-5, lOand 30 m lengths 

Additional Features 

• Common. field'Upgradeable 
software for aIIOEMV family receivers 
with OEM4 compatible commands 
and logs 

• Auxiliary strobe signals, including a 
configurable PPS output for time 
synchronization and mark inputs 

• Outputs to drive exterrul LEDs 


1 Typicai values l\iffc(nMrKe jpediicailons 
«ibjec( lo GPS system ch4i<Kte(Ktlc^ US 000 
operational ilegiadatioa lonospheik arvJ 
bopospherk coodilons. satellite eeometiy, 
baselne length, mulilpath effects end the 
presence of IntenOonal or urumenilonai 
aiterfeierKe sources. 

2 GPSonly. 

a expected eccurecy after static conwcrgeiKe. 

4 Slower data rales are expected for API cusiorrvers. 
The maxiinum data rale is dependent on the size 
of theappicailon 

5 Typlcaivalue. Noalmarucarephemertdesandno 
approximate podilon or time. 

6 Typkai value. Almanac and recent ephemef kies 
saved and apprcnimaie position and lime 
entered. 

7 Time accwacydoesnoiindudeblasesduetoRF 
or antenna deby. 

8 txpottlkendngresIricIsoperaUonloamaxImum 
of S14 metres per second. 

9 External CAN lianscelvef arvl user appHcatlon 
software required. 



1 SOCVNOVATCL (U.S. & Canada) or -403 29S-4900 
Europe ^4 (0) 1993 8S2-436 
SE Asia & Australia ^61 (0) 400 833 601 
saiestd novatei.com 

novatel.com 





Version 4 - SpedOcatlorvc subject to change without notice. 

C 7008 NuvA tel IrK. AH tights reserved. Printed In Canada. D095SS 
NovAteCand OEMV*are registered trademarks of NovAiel Inc. 
AtfVance*. GLIDE*.andRT-TO"*are tiadematksof NovAtelInc 
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4. Analog Devices ADIS16405 


□ ANALOG 
DEVICES 


Tri-Axis Inertial Sensor 
with Magnetometer 


ADIS16405 


FEATURES 

Tri-axis, digital gyroscope with digital range scaling 
±75Vsec, i 1 SOVsee, ±300Vsec settings 
Tri-axis, ± 18 g digital accelerometer 
Tri-axis, ±2.5 gauss digital magnetometer 
220 ms start-up time 

Factory-calibrated sensitivity, bias, and axial alignment 
Calibration temperature range: -40^C to -f 85^C 
Digitally controlled bias calibration 
Digitally controlled sample rate, up to 819.2 SPS 

External clock input enables sample rates up to 1200 SPS 
Digitally controlled filtering 
Programmable condition monitoring 
Auxiliary digital input/output 
Digitally activated self-test 
Programmable power management 
Embedded temperature sensor 
SPI-compatible serial interface 
Auxiliary, 12-bit ADC input and DAC output 
Single-supply operation: 4.75 V to 5.25 V 
2000 g shock survivability 
Operating temperature range: -40*C to±105*C 

APPLICATIONS 

Unmanned aerial vehicles 
Platform control 
Digital compassing 
Navigation 

GENERAL DESCRIPTION 

The ADISI6405 (Sensor*isacom{:4e(e inertial system that includes 
a tri-axis gyroscope, a tri-axis accelerometer, and a tri-axis mag¬ 
netometer. The ADIS16105 combines industry leading iMEMS* 
technology with signal conditioning that optimizes dynamic 
performance. The factory calibration characterizes each sensor 
for sensi! ivity, bias, alignment, and linear acceleration (gyroscope 
bias). As a result, each sensor has its own dynamic compensation 
for correct ion formulas that prov’kle accurate sensor measurements 
over a temperature range of -40®C to +85®C Tlie magnetometeis 
employ a self correction function to provide accurate bias 
performance over temperature as well. 

The ADISI6405 provides a simple, cost-effective method for 
integrating accurate, multi-axis inertial sensing into industrial 
systems, especially when compared with the complexity and 
investment associated with discrete designs. 


FUNCTIONAL BLOCK DIAGRAM 



All necessary motion testing and calibration are part of the produc¬ 
tion process at the factory, greatly reducing system integration time. 
Tight orthogonal alignment simplifies inertial frame alignment 
in navigation systems. An improved SPl interface and register 
structure provide faster data collection and contiguralion control. 
By using a compatible pinout and the same package as the 
ADlS1635xand ADIS1636X Emilies, upgrading to the ADIS16405 
requires only firmv^are changes to accommodate addit tonal sensors 
and register map i^datcs. 

This compact module is approximately 23 mm x 23 mm x 23 mm 
and provides a flexible connector interhkce, which enables mult^e 
mounting orientation options. 


R«v. 0 

w httg»*r ft nh«tiinUt l lmii ro wi nn 

cfthiRl paiMt tha mfl|r mU tram tu UM Sfwcffetfcn* todwgi miMmU nett* No 


On* T*chnoi«>gy W*y, P.O. Box 9106. Norwood, MA 02062-9106, U.S A. 
T«l: 781.3294700 www.analog.com 

Fax: 781.461.3113 C2009 AndegD*vic«s, Inc All rights r«s*rv«d. 
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ADIS16405 


SPECIFICATIONS 

Ta = -40*C to ♦85*C, VCC = 5.0 V, angular rale = 0*/sec, dynamic range = t300*/sec, 1 1 g, unless otherwise noted. 


Tabic 1. 


ParamotAT 

Test Conditions 

Min Typ Max 

Unit 

GYROSCOPtS 




Dynamic Range 


±300 ±350 

“/sec 

Initial Sensitivity 

Dynamic range s ±300Vsec 

0.0495 0.05 0.0505 

Vsec/LSB 


Dynamic range * ±1 SOVsec 

0.025 

‘/sec/LSB 


Dynamic range * ±75*/sec 

0.0125 

“/sec/LSB 

Sensitivity Temperature Coefficient 

-40"CsT*s+85*C 

±40 

ppm/*C 

Misalignment 

Axis-to-axis, & s 90° ideal 

±0.05 

Degrees 


Axis-to-firame (package) 

±0.5 

Degrees 

Nonlinearity 

Best f(t straight line 

0.1 

%ofFS 

Initial Bias Error 

1 a 

±3 

“/sec 

In-Run Bias Stability 

1o,SMPL_PRD = 0x01 

0.007 

“/sec 

Angular Random Walk 

lo.SMPL PRD = 0x01 

2.0 

“/Vhr 

Bias Temperature Coefficient 

-40“CsT*S+85“C 

±0.01 

“/sec/“C 

Linear Acceleration Effect on Bias 

Any axis, 1 o (MSC _aRL, Bit 7=1) 

0.05 

“/sec/g 

Bias Voltage Sensitivity 

VCC = 4.75Vto 5J5V 

0.32 

7sec/V 

Output Noise 

±300Vsec rartge, no filtering 

0.9 

•/sec rms 

Rate Noise Density 

f = 25Hi. ±300Vsec, no filtering 

0.05 

•/sec/VHz rms 

3 dB Bandwidth 


330 

Hz 

ACCELEROMETERS 




Dynamic Range 


±18 

9 

Initial Sensitivity 


3.285 3.33 3.38 

mg/LSB 

Sensitivity Temperature Coefficient 

-4(rCsT.s+85°C 

±50 

ppm/“C 

Misalignment 

Axis-to-axis, 0 = 90* Ideal 

OJ 

Degrees 


Axi$-to-firame (package) 

±0.5 

Degrees 

Nonlinearity 

Best fit straight line. ±17 g 

0.1 

%ofFS 

Initial Bias Error 

1 0 

±50 

mg 

In-Run Bias StaUlity 

1 0 

0.2 

mg 

Velocity Random Walk 

1 o 

0.2 

m/sec/Vhr 

Bias Temperature Coefficient 

-40“CsT*S+85*C 

±0.3 

mg/“C 

Bias Voltage Sensitivity 

VCC = 4.75Vto 5.25 V 

2.5 

mg/V 

Output Noise 

No filtering 

9 

mg rms 

Noise Density 

Noflttehng 

0.5 

mg/VHz rms 

3 dB Bandwidth 


330 

Hz 

MAGNETOMETER 




Dynamic Range 


±2.5 ±3.5 

gauss 

Initial Sensitivity 

2S“C 

0.49 0.5 O.Sl 

mgauss/LSB 

Sensitivity Temperature Coefficient 

25-^, 1 o 

600 


Axis Nonorthogonailty 

2 5°C, axis-to-axis 

0.25 


Axis Misalignment 

25*C axis-to-base plate and guide pins 

0.5 


Nonlinearity 

Best fit straight line 

0.5 

%ofF5 

Initial Bias Error 

2S‘^, 0 gauss stimulus 

±4 

mgauss 

Bias Temperature Coefficient 


0.5 

mgauss/“C 

Output Noise 

2 5*C no filtering 

1.25 

mgauss rms 

Ncrise Density 

2S“C, no filtering, rms 

0.066 

mgauss/VHz 

3 dB Bandwidth 


1540 

Hz 

TEMPERATURE SENSOR 




Scale Factor 

25"^. output = 0x0000 

0.14 

“C/LSB 
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ADIS16405 


Paramotor 

Test Conditions 

Min Typ Max 

Unit 

ADC INPUT 




Resolution 


12 

Bits 

Integral Nonlinearity 


±2 

LSB 

Differential Nonlinearity 


±1 

LSB 

Offset Error 


±4 

LSB 

Gain Error 


±2 

LSB 

Input Range 


0 3.3 

V 

Input Capacitance 

During acquisition 

20 

pF 

DAC OUTPUT 




Resolution 


12 

Bits 

Relative Accuracy 

Code 101 to Code4095,5 kO/lOO pF to GND 

±4 

LSB 

Differential Nonlinearity 


±1 

LSB 

Offset Error 


*5 

mV 

Gain Error 


tO.5 

% 

Output Range 


0 3.3 

V 

Output Impedance 


2 

D 

Output Settling Time 

SkO/lOOpF to GND 

10 

MS 

LOGIC INPUTS' 




Input High Voltage. Viw 


2.0 

V 

Input Low Voltage, Vi*«. 


0.8 

V 


CS signal to wake up from sleep mode 

0.55 

V 

« Wake-Up Pulse Width 


20 

MS 

Logic 11nput Current Imh 

V«=3.3V 

±0.2 ±10 

MA 

Logic 0 Input Current Imi 

< 

II 

o 

< 



All Pins Except R5T 


-40 -60 

ma 

^Pin 


-1 

mA 

Input Capacitance. Ch 


10 

pF 

DIGITAL OUTPUTS’ 




Output High Voltage, Voh 

isounas 1.6 mA 

2.4 

V 

Output Low Voltage, Va 

l&MKS 1.6 mA 

0.4 

V 

FLASH MEMORY 

Endurance^ 

10,000 

Cycles 

Data Retention* 

Tj = 85“C 

10 

Years 

FUNCTIONAL TIMES* 

Time until data is available 



Power-On Start-Up Time 

Normal mode, SMPL.PRD s 0x09 

220 

ms 


Low power mode, SMPL.PRD 2 OxOA 

290 

ms 

Reset Recovery Time 

Normal mode. SMPL.PRD s 0x09 

100 

ms 


Low power mode, SMPL PRD 2 OxOA 

170 

ms 

Sleep Mode Recovery Time 

Normal mode, SMPL _PRD s 0x09 

4 

ms 


Low power nrode, SMPL PRD 2 OxOA 

IS 


Flash Memory Test Time 

Normal mode, SMPL PRD s 0x09 

17 

ms 


Low power nxKle, SMPL_PRD 2 OxOA 

90 

ms 

Automatic Self-Test Time 

SMPL_PRD*0x01 

12 

ms 

CONVERSION RATE 

SMPL_PRD s 0x01 to OxFF 

0.413 819.2 

SPS 

Clock Accuracy 


±3 

% 

Sync Input Clock 


1.2 

kHz 

POWER SUPPLY 




Operating Voltage Range, VCC 


4.75 5.0 5.25 

V 

Power Supply Current 

Low power mode at 2 5*C 

45 

mA 


Normal mode 3125*^ 

70 

mA 


Sleep mode at2SX 

600 

ma 


' The digitdl IA!> signals are driven by an inicvnal 3.3 V supply, and the inputs are $ V tderanl. 

^ Endurance isgualified as per lEOEC Starsdard 2Z Method A117, artd measured at -40*C +2S*C 48S*C and +12S*C 

' The data retention lilctime equivalent is at a function temperature (Tj) of 8SX as per 3EHC $ tarsdard 22. Method At 17. Data retention lifeOmc decreases with fimction 
tcmpcrative. 

* These times do not indude thermal setilirsg and intcrrral filter response limes (330 Hz bandwidth), which may affect ovcsali accuracy. 
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TIMING SPECIFICATIONS 

Ta = 25*C, VCC = 5 V, unless otherwise noted. 

Table 2. 


Parameter 

fsux 

t$Ull 

tRMAAtl 

t<S 

tOAV 

tosu 

tOFO 

tsauRf tsaie 
tof, toe 
tvs 



Normal Mode Low Power Mode 

(SMPL PROS 0x09) (SMPL PRD 2 0x0A) Burst Mode 


Description 

Min' Typ Max 

Min' Typ Max 

Min' Typ Max 

Unit 


0.01 2.0 

0.01 0.3 

0.01 1.0 

MHz 

Stall period between data 

9 

75 

1/fs<iK 


Read rate 

40 

150 


ps 

Chip select to clock edge 

4S.8 

48.8 

48.8 

rts 

DOUT valid after SCLK edge 

100 

100 

100 

ns 

DIN setup time before SCLK rising edge 

24.4 

24.4 

24.4 

ns 

DIN hold time after SCLK rising edge 

48.8 

48.8 

48.8 

ns 

SCLK rise/fall times 

5 123 

5 12.5 

5 123 

ns 

DOUT rise/fall times 

5 12.5 

5 123 

5 12.5 

ns 

CS high after SCLK edge 

5 

5 

5 

ns 

Input sync pulse width 

5 



MS 

Input sync to data ready output 

600 



MS 

Input sync period 

833 



MS 


'Guaranteed b/ design and characterizatioa but not tested in production. 





Figured Input Clock Timing Diagram 
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ABSOLUTE MAXIMUM RATINGS 


Tables. 


Parameter 

Rating 

Acceleration 


Any Axi 5 » Unpowered 

2000 9 

Any Axis, Powered 

2 OOO 9 

VCCtoGND 

-0.3 V to +6.0 V 

Digital Input Voltage to GND 

-a3Vto+53V 

Digital Output Voltage to GND 

-0.3VtoVCC + 0.3V 

Analog Input to GND 

-0.3 V to+3.6 V 

Operating Temperature Range 

-40’C to+105"C 

Storage Temperature Range 

-65"Cto+125’C'' 


' Extended exposure to temperatures outside the specified remperaiixe 
range of -40*C to 4^10S*Ccan adversely affect the accuracy of the factory 
calibration. For best accuracy, store the parts within the specified operating 
range of -40*C to r-10S*C 

^ Although the device is capable of withstanding short-term eigKKure to 
1S0*C long-term exposure threatens internal mechanical integrity. 


Stresses above those listed under Absolute Maximum Rat ings 
may cause permanent damage to the device. This is a stress 
rating only; functional operation of the device at these or any 
other conditions above those indicated in the operational 
section of this specification is not implied Exposure to absolute 
maximum rating conditions for extended periods may affect 
dcWcc reliability. 


Table 4 . Package Characteristics 


Package Type 

6ja 

6jc 

Device Weight 

24-Lead Module 

39.8XAV 1 

14.2"C/W 

16 grams 


ESD CAUTION 




ESO (■toctiusUlk diictuirg*) Mnsitiv* d«v<c«. 

Charg^ devices end circuit boards can discharge 
without detection. Although this produa features 
patented or proprietary protection circuitry, damage 
nnay oaur on dmdces subjected to high energy ESO. 
There fo r e , proper ESO precautions should be taken to 
avoid performance degradation or loss of functionality. 
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PIN CONFIGURATION AND FUNCTION DESCRIPTIONS 


ADIS16405 
TOP VIEW 
(Not to Scale) 
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NOTES 

1. THIS VIEW REPRESENTS THE TOP VIEW OF THE MATING CONNECTOR. 

2. WHEN CONNECTED TO THE ADIS164t)5, THE PINS WILL NOT BE VISIBLE. ^ 

3. MATING CONNECTOR: SAMTEC CLM-112-02 OR EQUIVALENT. Z 

4. DNC = DO NOT CONNECT. | 

Figures. PinConftguration 



Figure 6. Axial Orientoffon 


Table 5. Pin Function Descri 

>tions 

Pin No. 

Mnemonic 

Type’ 

Description 

1 

DI03 

I/O 

Configurable Digital Input/Output. 

2 

DI04/CLKIN 

I/O 

Configurable Digital Input/Output or Sync Clock Input. 

16,17,18,19, 22, 23, 24 

DNC 

N/A 

Do Not Connect. 

3 

SCLK 

1 

SPI Serial Clock. 

4 

DOUT 

0 

SPI Data Output. Clocks output on SCLK falling edge. 

5 

DIN 

1 

SPI Data Input. Clocks input on SCLK rising edge. 

6 

CS 

1 

SPI Chip Select. 

7 

DIOl 

I/O 

Configurable Digital Input/Output. 

8 

RST 

1 

Reset. 

9 

DI02 

I/O 

Configurable Digital Input/Output. 

10,11,12 

VCC 

s 

Power Supply. 

13.14.15 

GND 

s 

Power Ground. 

20 

AUX_DAC 

0 

Auxiliary, 12-Bit DAC Output. 

21 

AUX ADC 

1 

Auxiliary, 12-Bit ADC input. 


' S is supply, O is output, I is input, N/A is not applicable. 
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ADIS16405 


THEORY OF OPERATION 

BASIC OPERATION 

The ADISI6405 is an autonomous sensor system that starts up 
after a valid power supply voltage is applied and then begins 
producing inertial measurement data at the factory default sample 
rate of 819^ SPS. After each sample cycle, the sensor data loads 
into the output registers and DIO 1 pulses, providing a new data 
ready control signal for driving system -level internet service 
routines. In a typical system, a master processor accesses the 
output data registers throu^ the SPl internee, using the hook 
up shown in Figure 9. Table 6 provides a generic functional 
description for each pin on the master processor. Table 7 
describes the typical master processor settings normally found 
in a configuration register and used for communicating with 
thcADIS16405. 


UO UNE8 ARE COMPATIBLE VMTH 
3 3V OR 6V LOGIC LEVELS 



'rable6. Generic .Master Processor Pin Names and Functions 


Pin Name 

Function 

SS 

Slave select 

IRQ 

Interrupt request 

MOSI 

Master output, slave input 

MISO 

Master input, slave output 

SCLK 

Serial clock 


Table 7. Generic Master Processor SPl Settings 


Processor Setting 

Description 

Master 

The ADISI 6405 operates as a slave. 

SCLK Rate £2 MHz' 

Normal mode. SMPL_PRDt7:0] s 0x08. 

CPOL=l 

Clock polarity. 

CPHA=1 

Clock phase. 

MSB-First 

Bit sequence. 

16-Bit 

Shift register/data length. 


' For butsi mode SO.Kraie s 1 MKc For tow power mode SCLK rales 300 Idti. 


The user registers provide addressing for all inpuf/output 
operations on the SPl interface. Each 16-bit register has two 
7-bit addresses: one for its i^^per byte and one for its lower byte. 


Table 8 lists the lower byte address for each register, and F^re 10 
shows the generic bit assignments. 

|.s|l4[.3|»|..|..|.|.|r|.|.|.[3|2[.|.| 

I -UPPBIBYTB-1 LOWBI BYTE-1 | 

Figufe W. Output Register Bit Assignments 

READING SENSOR DATA 

Although the ADISI6405 produces data independently, it oper¬ 
ates as an SPl slave device that communicates with system (master) 
processors using the 16-bit segments displayed in Figvue 11. 
Individual register reads require two such 16-bit sequences. The 
first 16 bit sequence provides the read command bit (RAV = 0) 
and the target register address (A6 to AO). The second sequence 
transmits the register contents (D15 to DO) on the DOUT line. 
For example, if DIN = OxOAOO, then the content of XACCL_OUT 
shifts out on the DOUT line during the next 16-bil sequence. 

Tlie SPl operates in lull duplex mode, whidi means that the master 
processor can read the output data from DOUT while using the 
same SCLK pulses to transmit the next target address on DIN. 

DEVICE CONFIGURATION 

The user register memory map (Table 8) identifies configuration 
registers with either a W or R/W. Configuration commands also 
use the bit sequence displayed in Figure 12. If the MSB is equal to i, 
the last ei^t bits (DC7 to DCO) in the DIN sequence load into 
the memory address associated widi die address bits (AS to AO). 
For example, if DIN s OxAl IF, then 0x1 F loads into Address 
Location 0x21 (XACCL.OFF, upper byte) at the conclusion of 
the dat a frame. 

Most of the registers have a backup location in nonvolatile flash 
memory. The master processor must manage the backup function. 
Set GLOB_CMD[3l = 1 (DIN = 0xBE04) to execute a manual 
tladi update (backup) operation, which copies the user registers 
into their respective flash memory locations. This operation 
takes SO ms and requires the power supply voltage to be within 
the specitied limit to complete properly. The FLASH..CNT register 
provides a running count of these events tor managing the long 
term reliability of the flash memory. 

BURST MODE DATA COLLECTION 

Blits! mode data collection offers a more efficient method tor 
collecting data from the AD1SI640S. In sequential data cycles 
(each separated by one SCLK period), all output registers clock 
out on DOUT. This sequence starts when die DIN sequence is 
9011 1110 00(K) 0000 (OxSEOO). Ncja, the contents of each output 
register are output from DOUT, starting with SUPPLY_OUT 
and ending w'hh AUX_ADC (see Figure 12). The addressing 
sequence shown in Table 8 determines the order of the outputs 
in burst mode. 
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Tabic 8. User Register Mcmor>' Map 


Name 


Flash Backup 


Default 

Function 

Bit Assignments 

FLASH.CNT 

R 

Yes 

QkOO 

N/A 

Flash memory write count 

N/A 

SUPPLY.OUT 

R 

No 

0x02 

N/A 

Power supply measurement 

Table 9 

XGYRO.OUT 

R 

No 

0x04 

N/A 

X-axis gyroscope output 

Table 9 

YGYRO_OUT 

R 

No 

0x06 

N/A 

Y-axis gyroscope output 

Table 9 

ZGYRO.OUT 

R 

No 

0x08 

N/A 

Z-axis gyroscope output 

Table 9 

XACCL^OUT 

R 

No 

OxOA 

N/A 

X-axis accelerometer output 

Table 9 

YACCL^OUT 

R 

No 

OxOC 

N/A 

Y-axis accelerometer output 

Table 9 

2ACCL OUT 

R 

No 

OxOE 

N/A 

Z-axis accelerometer output 

Tables 

XMAGNOUT 

R 

No 

0x10 

N/A 

X-axis magnetometer measurement 

Tables 

YMAGN OUT 

R 

No 

0x12 

N/A 

Y-axis magnetometer measurement 

Tables 

ZMAGN OUT 

R 

No 

0x14 

N/A 

Z-axis magnetometer measurement 

Tables 

TEMP_OUT 

R 

No 

0x16 

N/A 

Temperature output 

Tables 

AUXADC 

R 

No 

0x18 

N/A 

Auxiliary ADC measurement 

Tables 

XGYRO_OFF 

mi 

Yes 

OxIA 

0x0000 

X-axis gyroscope bias offset factor 

Table 10 

YGYRO.OFF 

RAW 

Yes 

OxlC 

0x0000 

Y-axis gyroscope bias offset factor 

Table 10 

ZGYRO_OFF 

R/W 

Yes 

OxlE 

0x0000 

Z-axis gyroscope bias offset factor 

Table 10 

XACCL.OFF 

R/W 

Yes 

0x20 

0x0000 

X-axis acceleration bias offset factor 

Table 11 

YACCL.OFF 

R/W 

Yes 

0x22 

0x0000 

Y-axis acceleration bias offset factor 

Table 11 

ZACCL.OFF 

RAW 

Yes 

0x24 

0x0000 

Z-axis acceleration bias offset factor 

Table 11 

XMAGN HIF 

RAW 

Yes 

0x26 

0x0000 

X-axis magnetometer, hard-iron factor 

Table 12 

YMAGN.HIF 

RAW 

Yes 

0x28 

0x0000 

Y-axis magnetometer, hard-iron factor 

Table 12 

ZMAGN.HIF 

RAW 

Yes 

Qx2A 

0x0000 

Z-axis magnetometer, hard-iron factor 

Table 12 

XMAGN_SIF 

RAW 

Yes 

Gx2C 

0x0000 

X-axis magnetometer, soft-iron factor 

Table 13 

YMAGN.SIF 

RAW 

Yes 

0x2E 

0x0000 

Y-axis magnetometer, soft-iron factor 

Table 13 

ZMAGN_SIF 

RAW 

Yes 

0x30 

0x0000 

Z-axis magnetometer, soft-iron factor 

Table 13 

GPIO.CTRL 

RAW 

No 

0x32 

0x0000 

Auxiliary digital input/output control 

Table 18 

MSC.CTRL 

RAW 

Yes 

0x34 

0x0006 

Miscellaneous control 

Table IS 

SMPL. PRD 

RAW 

Yes 

0x36 

0x0001 

Internal sample period (rate) control 

Table 15 

SENS AVG 

RAW 

Yes 

0x38 

0x0402 

Dynamic range and digital filter control 

Table 17 

SLP_CNT 

W 

No 

Qx3A 

0x0000 

Sleep mode control 

Table 16 

DIAG.STAT 

R 

No 

0x3C 

0x0000 

System status 

Table 23 

GLOB.CMD 

W 

N/A 

0x3E 

0x0000 

System command 

Table 14 

ALM.MAGl 

RAW 

Yes 

0x40 

0x0000 

Alarm 1 amplitude threshold 

Table 25 

ALM_MAG2 

RAW 

Yes 

0x42 

0x0000 

Alarm 2 amplitude threshold 

Table 25 

ALM_SMPL1 

RAW 

Yes 

0x44 

QxOOOO 

Alarm 1 sample size 

Table 26 

ALM_SMPL2 

RAW 

Yes 

0x46 

0x0000 

Alarm 2 sample size 

Table 26 

ALM_CTRL 

RAW 

Yes 

0x48 

0x0000 

Alarm corttrol 

Table 24 

AUX_DAC 

RAW 

No 

0x4A 

0x0000 

Auxiliary DAC data 

Table 20 


' Each re 9 lsicr contains two bytes. The address of the Iowa byte is displayed. The address of the upper byte is equal to the address of the lower byte plus T. 


=*“1_rti_ 

DIN (^jf Ae Has ITm I ai | ao Ioc 7 loce |dc6 |dc4 |dc3 loca loci |dco A6 I Ae 

DOUT —^^^01Sj^14jJouJ^D12j[01lJ^D10jLDOjL08jQ^j[_DeJ_06jL04j_raJ_ra_^01_^00 }-k: ^OWjDUj^Dt^ 
NOTE$ 

1. OOUT BITS ARE BASED ON THE PREVIOUS ie«IT SEQUENCE |R> 0). 

Figure 1 lOutput Regifter Bit Assignmerra 
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OUTPUT DATA REGISTERS 

Figure 6 provides ihc posilivc measurement direction for eacli 
gyroscope, accelerometer, and magnetometer. Table 9 provides 
the condguration and scale factor for each output data register 
in the ADIS16405. All inertial sensor outputs are Mbitsinler^th 
and are in twos complement format, which means that 0x0000 
is equal to 0 LSB, 0x0001 is equal to +1 LSB, and OxBFFF is equal 
to - I LSB. The following is an example of how to calculate the 
sensor measurement from the XGYRO_OUT: 

XGVRO.Ot/r =0x3B4A 

0x000 - 0X3B-1A = -0x04B6 = (4 x 256 +11x 16 +6)- 
0x04B6--1206 LSB 



Figure 13. Ckapuf Reciter Bit Aiugnmenis 


Auxiliary ADC 


The AUX_ ADC register provides access to the auxiliary ADC 
input channel. The ADC is a 12 -bit successive approximation 
converter that has an equivalent ii^ut circuit to the one shown 
in Figure 14. The maximum input is 3.3 V. Tltc ESD protection 
diodes can handle 10 mA without causing irreversible damage. 
The on resistance (Rl)of the switch has a typical value of 100 D. 
The sampling capacitor, C2, has a typical value of 16 pF. 
vcc 


Rate = 0.05*/sec x (-1206) = -603*/sec 
Therefore, an .XGYRO_OUT output of0x3B4A corresponds to 
a clockwise rotation about the z-axis (see Figure 6) of 60.3*/sec 
when looking at the top of the package. 

Tabic 9. Output Data Register Formats 


Register 

Bit, 

Format 

Scale 

SUPPLY.OUT 

14 

Binary, 5 V = 0x0814 

242 mV 

XGYRO.OUT’ 

14 

Twos complement 

O.OSVsec 

YGYRO.OUT' 

14 

Twos complement 

O.OSVsec 

ZGYRO_OUT' 

14 

Twos complement 

0.05‘’/sec 

XACCL.OLTT 

14 

Twos complement 

lOmj 

YACCL.OLIT 

14 

Twos complement 

lOm^ 

ZACCL.OUT 

14 

Twos complen'>ent 

10 mg 

XMAGN_OUT 

14 

Twos complement 

0.5 mgauss 

YMAGN^OUT 

14 

Twos complement 

0.5 mgauss 

ZMAGN_OUT 

14 

Twos complement 

0-5 mgauss 

TEMP.OUT' 

12 

Twos complement 

0.14‘C 

AUX_ADC 

12 

Binary,! V = 0x04D9 

0.81 mV 


' AsumM that the «cjting K set to ±300*/sec This betor scales with the ranoe. 
^ The typical output for this register at is OkOOOO. 


Each output data register uses the bit assignments shown in 
Figure 13. The ND flag indicates that unread data resides in the 
output data registers. Tliis Hag clears and retiu+is to 0 during an 
output register read sequence. It returns to 1 atterthc next internal 
sample iq>dates the registers with new data. The EA llag indicates 
that one of the error 6i^s in the DIAG. STAT register (see TaUe 23) 
is active (true). The remaining 14 bits are for data. 



Figure 14. Bquivaienl Analog Input Circuit 
(Corwersion Phase: Switch Open. Track Phase: Switch Closed) 


CALIBRATION 


Manual Bias Co//draf/on 


The bias offset registers in Table 10. Table 11, and Table 12 
(hard-iron correction for magnetometer) provide a manual 
adjustment function for the output of each sensor. For example, 
if XGYRO.OFF equals 0xlFF6. the XGYRO.OUT offset shifts 
by -10 LSB, or -0.125*/s<c. The DIN command for the upper 
byte is DIN = 0x9B IF; for the lower byte, DIN = 0x9AF6. 


Table 10. XGYRO_OFF, YGYRO^OFF, ZGYRO.OFF 


Bits 

Desmption 


Not used. 

[12«1 

Data bits. Twos complement. 0.01257$ec per LSB. 
Typical adjustment range » ±507sec. 


Table II . XACCL. OFF, YACCL. OFF. ZACCL .OFF 


Bits 

Description 

[15:12) 

Not used. 

[110] 

Data bits. Twos complement, 3-3 mg/LSB. 

Typical adjustment range = ±6.75 g. 


Table 12. XMAGNJIIF, YMAGN_MIF,ZMAGN_H1F 


Bits 

Description 

115:14) 

Not used. 

[13*) 

Data bits. Twos complement, 0.5 mgauss/LSB. 

Typical adjustn^nt range = ±4 gauss. 
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Magnetometer Soft-Iron Correction (Scale Factor) 

The sott iron correction factor for the magnetometer provides 
opportunity to citangc the scale (actor for eadi individual axis. 

Table 13. XMAGN SIF, YMAGN,S1F,ZMAGN SIF 


Bits 

Description 


Not used. 


Data bits. Binary, linear scale adjustment factor 
between OxOOCA (Ox) and OxSFFF (2x). 


Gyroscope Automatic Bias Null Calibration 

Set GLOB_CMD[0]= 1 (DIN = 0xBE01)toexeculelhislimction, 
which measures the gyroscope outputs and then loads the gyro¬ 
scope offset registers with the opposite values to provide a quick 
bias calibration. Tlien, all sensor data resets to 0, and tlie flash 
memory updates automatically within 50 ms(sceTablc 14). 
Gyroscope Precision Automatic Bias Null Calibration 
Set GLOB_CMD(4l = I {DIN = OxBEI 0) to execute this liinction, 
which lakes the sensor offline for 30 sec while it collects a set of 
gyroscope data and calculates a more accurate bias correction 
factor for each gyroscope. Once calculated^ the correct ion factor 
loads into the three gyroscope offset registers, all sensor data 
resets to 0. and the flash memory updates automatically within 
50 ms (sec Table 14). 

Restoring Factory Calibration 

Set GLOB_CMDl 1) = I (DIN = OxBE02) to execute this timetion, 
which resets each user calibration register (sec Table 10, Table 11, 
and Figure 12) to 0x0000, resets all sensor data to 0, and auto- 
maticaUy updates the flash memory within 50 ms (see Table 14). 
Linear Acceleration Bias Compensation (Gyroscope) 

Set MSC_Cn'RL[7| = I (DIN = 0xB486) to enable correction for 
tow frequency acceleration influences on gyroscope bias. Note 
that the DIN sequence also preserves the factory default condi¬ 
tion for the data ready function (see Table 19). 

OPERATIONAL CONTROL 
Global Commands 

The GLOB_CMD register provides trigger bits for several useful 
functions. Setting tlie assigned bit to 1 starts eadt operation, and the 
bit returns to 0 after conviction. For example, set GLOB_CMD(7] 
* 1 (DIN * OxBESO) to execute a software reset, which stops the 
sensor operation and runs the device throu^ its start-up sequence. 
This includes loading the control registers with their respect we 
flash memory locations prior to producing new data. Reading 
the GLOB_CMD registers (DIN = 0x3E00) startsihe burst 
mode read sequence. 


Table 14. GLOB_CMD 


Bin 

Description 

[15:8] 

Not used 

[7] 

Software reset command 

[6:51 

Not used 

[41 

Precision autonull command 

[31 

Flash updatecommand 

[21 

Auxiliary DAC data latch 

[11 

Factory calibration restore command 

[01 

Autonull command 


Internal Sample Rate 

The ADIS16405 performs best when the sample rate is set to the 
factory default setting of SI 9.2 SPS. For applications that require 
lower sample rates, the S.VIPL PRD register controls the AD1S16405 
internal sample (sec Table 15), and the following rclationsliip 
produces the sample rate: 


ts 

= h X Nj + 1 

Table IS. SMPL PRD 

Bits 

Description 

[15:81 

Not used 

m 

Time base (ts) 

0 = 0.61035 ms,l =ia921 ms 

[601 

Increment setting (N;) 

internal sample period »ts ~ te x N^ + 1 


For example, set SMPL_PRD(7K)] ^OxOA (DIN = 0xB60A) for 
an internal sample period of 6.7 ms and a sample rale of 149 SPS. 

Power Management 

Setting SMPL. PRD 2 OxOA also sets the sensor in low power 
mode. For systems that require lower pKiwer dissipation, in system 
characterization helps users quality the associated pHrrformance 
trade-offs. In addition to sensor performance, this mode affects 
SPI data rates (see Table 2). Two sleep mode options are listed 
in Table 16. Set SLP_CNTI8] = I (DIN = OxBBOl) to start the 
indefinite sleep mode, which requires a CS assertion (high to 
low), reset, or p>owercyde to wake up. Set SLP„CNT[7K)J = 0x6>l 
(DIN - OxBAM) to pul the ADIS16405 to sleep for 50 sec, as 
an example of the programmable sleep time option. 

Tabic 16. $IP_CNT _ 

Bits Description _ 

(15:9] Not used 

[S] Indefinite sleep mode, set to 1 

[7K)] Programmable sleep tinne bits, 0.5 sec/LSB 


Digital Filtering 

Programmable low pass filtering provides additional Ofvoriunity 
for noise reduction on the inertial sensor outputs. This f ilter 
contains two cascaded averaging filters that provide a Bartlett 
window, FIR filter response (see Figure 15). SENS_AVG[2:0l 
controb the number of taps in each averaging stage. For example. 
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SENS_AVG(2K)] = 110 sets each stage tap to 64. The total 
number of taps in the filter is equal to 2N i 1. 



0.001 0 J )1 0.1 1 I 

FREQUENCY (Tf,) | 

figure 15. Btulhrti Window FK Frequency Response 
(Phase = N Samples^ 

Dynamic Range 

There are three dynamic range settings for the gyroscope; ±75*/scc» 
±l50*/sec, and iSOO^/scc. The lower dynamic range settings 
(±75*/sec and ± 1507sec) limit the minimum filter tap sizes to 
maintain the resolution as the measurement range decreases. 
The recommended order for programming the SENS.AVG 
register is upper byte (sensitivity), followed by lower byte 
(filtering). For example, set SENS_AVGlI0:8l = 010 (DIN 
0xB9O2) for a measurement range of ±1507sec, and then set 
SENS_AVG(2K)I = 110 ( DIN = 0xB806) for 64 taps per stage 
(129 taps overall). 


Tabic 17. SENS.AVG 


Bit. 

Settinfp 

Desoiption 

[15:11] 


Not used 

110:8) 


AAeasurement range (sensitivity) selection 


100 

l300”/sec (default condition) 


010 

±1507sec, fitter taps 2 4 (Bits[2:0] 2 0x02) 


001 

±757sec, filter taps 2 16(Brts[2;0] 2 0x04) 

[7:31 


Not used 

[2:0) 


Number of taps in each stage N » 2*^ 


INPUT/OUTPUT FUNCTIONS 

General-Purpose I/O 

DIOl. D102. DI03, and DI04 are configurable, general 
purpose I/O lines lliat serve muitqHe purposes according to the 
following control register priority: MSC.CTRU ALM_CTRL, 
and GPIO.CTRL. For example, set GPIO_CTRL = OxOSOC 
(DIN = OxBSOS. and then 0xB40C) to set DIOl and DI02 as 
inputs and DI03 and DI04 as outputs, with DIOS set low and 
DI04 set high. 


Tabic 18. GPIO.CTRL 


Bits 

Desoription 

[15:121 

Not used 

[111 

General-Purpose I/O Line4 {DI04) data level 

[101 

General-Purpose I/O Line 3 (DI03) data level 

[91 

General-Purpose I/O Line 2 (DI02) data level 

[81 

General-Purpose I/O Line 1 (DIOl) data level 

[7S11 

Not used 

[3] 

General-Purpose I/O Line 4 (DI04). direction control 

1 = output. 0 = input 

[2] 

General-Purpose I/O Line 3 (DI03). direction control 

1 s output. Os input 

[11 

General-Purpose I/O Line 2 (DI02). direction control 

1 3 output, 0 s input 

[0! 

General-Purpose I/O Line 1 (DIOl). direction control 

1 3 output, 0 = input 


Input Clock Configuration 

The input clock allows for external control over sampling in the 
ADIS16405. Set GP10_CTRL(3l = 0 (DIN = 0x0B200) and 
SMPL_PRDl7;0l = 0x00 (DIN - 0x8600) to enable this fonction. 
See Table 2 and Figure 4 for timing information. 

Data Ready I/O Indicator 

The factory default sets DIOl as a positive data ready indicator 
signal. The MSC_CTRL(2.'0] register provides configuration 
options for changing this. For example, set MSC,CTRL(2K)1 = 
100 (DIN = 0x8404) to change the polarity of the data ready 
signal for interrupt inputs that require negative logic inputs for 
activation. The pulse width will be between 100 ps and 200 ps 
over all conditions. 


Tabic 19. MSC CTRL 


Bits 

Description 

[15:12] 

Not used 

[111 

Memory test (clears on completion) 

1 = enabled. 0 s disabled 

[10] 

Internal self-test enable (clears on completion) 

1 = enabled, 0 = disabled 

[91 

Manual self-test, negative stimulus 

1 s enabled, 0 s disabled 

[81 

Manual self-test positive stimulus 

1 s enabled. Os disabled 

[71 

Linear acceleration bias compensation for gyroscopes 

1 3 enabled, 0 » disabled 

[61 

Linear accelerometer origin alignment 

1 3 enabled, 0 3 disabled 

[5:31 

Not used 

[21 

Data ready enable 

1 3 enabled, 0 = disabled 

[11 

Data ready polarity 

1 3 active high, 0 3 aaive low 

[01 

Data ready line select 

1 =DIO2,0=DIOl 
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AuxlUary DAC 

The 12 -bit AUX_DAC line can drive its output to within 5 mV 
of the ground retercnce when it is not sinking current. As the 
output approadies 0 V, the linearity begins to degrade (-100 LSB 
beginning point). As the sink current increases, the nonlinear 
range increases. The DAC latch command moves the values of 
the AUX_DAC register into the DAC input register, enabling 
both bytes to take effect at the same time. 



DIAGNOSTICS 

Self-Test 

The self-test function offers the opportunity to verify the 
mechanical integrity of each MEMS sensor. It applies an 
electrostatic force to each sensor element, which results in 
mechanical displacement Uiat simulates a response to actual 
motion. Table I lists the e;q>ected response for eacli sensor, 
which provides pass/fail criteria. Set MSC_CTRL[I0] = 1 (DIN 
s 0xB5O4) to run the internal self-test routine, which exercises 
all inertia] sensors, measures each re^nse, makes pass/fail 
decisions, and reports them to error flags in the DIAG_STAT 
register. MSC_CTRL(10) resets itself to 0 after completing the 
routine. MSC_CrrRI,l9:8l (DIN = 0xB502 orOxBSOl) provides 
manual control over the self-test function. Table 22 shows an 
example test flow for using this option to check the x-axis 
gyroscop>e. Zero motion provides results that are more reliable. 
The settings in Table 22 arc flexible and provide opportunity for 
optimization around speed and noise influence. For example, 
using fewer filtering taps decreases delay times but increases the 
opportunity for noise influence. 

Memory Test 

Setting MSC_CTRL|l ll = I (DIN = OxBSOS) performs a 
checksum verification of the flash memory locations. The 
pass/fail result loads into the DlAG_STAT[6l register. 

Sfofus 

The error flags provide indicator functions for common system- 
level issues. All of the flags clear (set to 0) after each DIAG _STAT 
register read cycle. If an error condition remains, the error flag 
returns to 1 during the next sample cycle. DIAG.STATll.-Oj 
does not require a read of this register to return to 0. 


Table 22. Manual Self-Test Example Sequence 


DIN 

Description 

0XB601 

SMPL_PRDt7:0I = 0x01, sample rate = 819.2 SPS. 

0xB904 

SENS_AVG[15:8] = 0x04, gyroscope range = ±3(X)7sec. 

0xB802 

SENS .AVG[7K)] s 0x02, four-tap averaging filter. 

Delay s 50 ms. 

0x0400 

Read XGYRO.OUT. 

0XB502 

MSC,CTRU91 * 1. gyroscope negative self-test. 

Delay « 50 ms. 

0x0400 

Read XGYRO.OUT. 

Calculate the positive change from the first reading 
to the second reading of XGYRO_OUT, and check to 
make sure the change is within the positive self-test 
response range specified in Table 1. 

OxBSOI 

MSC. CTRL[9:d] sOI. gyroscope/accelerometer 
positive self-test. 

Delay = 50 ms. 

0x0400 

Read XGYRO.OUT. 

Calculate the negative change from the first reading 
to thethird reading of XGYRO_OLrr, and check to 
make sure the change is within the positive self-test 
response range specified in Table 1 

0x8500 

MSC_CTRLI15:8]=0x00. 


3. DIAG_STATBU Descriptions 
Description 

Z'dxis accelerometer $etf*te$t failure (1 = fail, 0 = pass) 
Y-axis accelerometer self-test failure (1 * fail 0 * pass) 
X-axis accelerometer selfrtest failure (1 « fail, 0 ~ pass) 
X-axis gyroscope self-test failure (1 s fail, 0 s pass) 
Y-axis gyroscope self-test failure (1 s fail, 0 s pass) 
Z-axis gyroscope self-test failure (1 s fail, 0 s pass) 
Alarm 2 status (1 = active, 0 = inaaive) 

Alarm 1 status (1 = active, 0 = inactive) 

Not used 

Flash test, checksum flag (1 s fail, 0 = pass) 

Seif-test diagnostic error flag (1 s fail, 0 = pass) 

Sensor ovetrange (1 = fail, 0 = pass) 

SPI communication failure (1 s fail, 0 s pass) 

Flash update failure (1 s fait 0 s pass) 

Power supply above S.25 V 

(1 s power supply 2 5.25 V, 0 « power supply s 5.25 V) 
Power supply below 4.75 V 

(1 = power supply i 4.75 V. 0 = power supply 2 4.75 V) 
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Alarm Registers 

The alarm function provides monitoring for two independent 
conditions. The ALM_CTRL register provides control inputs 
for data source, data liltering {prior to comparison), static 
comparison, dynamic rate-of-change comparison, and output 
indicator configurations. The ALM .MAGx regi^ers establish the 
trigger threshold and polarity configurations. 

Table 27 gives an example of how’ to configure a static alarm. 
The ALM_SMPLx registers provide the number of samples to 
use in the dynamic rate-of change configuration. The period 
equals the number in the ALM_SMPLx register multiplied by 
the sample period time, which is established by the SMPL_PRD 
register. See Table 28 for an example of how to configure the 
sensor for ihistype of function. 


Table 24. ALM CTRL Bit Designations 


Bits 

Sattings 

Dascription 

H5:12] 


Alarm 2 source selection 


0000 

Disable 


0001 

Power supply output 


0010 

X-axis gyroscope output 


0011 

Y-axis gyroscope output 


0100 

Z-axts gyroscope output 


0101 

X*axis accelerometer output 


0110 

Y>axi$ accelerometer output 


0111 

Z*axis accelerometer output 


1000 

X*axis magnetometer output 


1001 

Y-axis magnetometer output 


1010 

Z-axis magnetometer output 


1011 

Gyroscope temperature output 


1100 

Auxiliary ADC input 

tn:8J 


Alarm 1 source selection (same as Alarm 2) 

m 


ftate-of-char>ge (ROC) enable for Alarm 2 

1 = rate of change 0 = static level 

[6] 


Rate-of-char>ge (ROC) enable for Alarm 1 

1 s rate of change 0 s static level 

(51 


Not used 

(41 


Comparison data filter setting' 

1 s filtered data, 0 s unfiltered data 

(31 


Not used 

(21 


Alarm output enable 

1 s enabled, 0 s disabled 

(11 


Alarm output polarity 

1 s active high, 0 s active low 

(01 


Alarm output line select 
lsDtO2.0«D)Ol 


' Incline outputs always use tiltefed data in this comparison. 


Tabic 25. AI..M_MAG1, AIJVI_MAC:2 

Bits 

Description 

(151 

Comparison polarity 

1 = greater than. 0 = less than 

(141 

Not used 

(133)1 

Data bits that match the format of the trigger source 
selection 

Table 26. AL.M_SMPL1, AIAI_SMPL2 

Bits 

Descriptien 

[15:«1 

Not used 

(73)1 

Data bits: number of samples (both 0x(K) and 0x01 ~ 1) 


Table 27. Alarm Cx)nriguration Example 1 


DIN 

Description 

OxAFSS, 

ALM_CTRL = 0x5517. 

0xA£17 

Alarm 1 input = XACCL_OLrr. 

Alarm 2 input = XACCL_OUT. 

Static level comparison, filtered data. 

DI02 output indicator, positive polarity. 

0xA783, 

ALM MAGI =0x8341. 

0xA641 

Alarm 1 is true if XACCL.OUT > 0.5 g. 

0xA93C. 

ALM MAG2= 0x3CBF. 

OxASBF 

Alarm 2 is true if XACCL_OLrr < -0.5 g. 


Table 28. Alarm (Configuration Example 2 


DIN 

Description 

0xAF76^ 

0xA£87 

ALM_aRL = 0x7687. 

Alarm 1 input = ZACCL. OUT. 

Alarm 2 Input * YACCL.OUT. 

Rate of change comparison, unfiltered data. 

DI02 output indicator. F>ositive polarity. 

0xB601 

SMPL PRD = 0x0001. 

Sample rate = 819.2 SPS. 

OxABOS 

ALM SMPLl = 0x0008. 

Alarm 1 rate of change period = 9.77 ms. 

OxACSO 

ALM_SMPL2= 0x0050. 

Alarm 2 rate of change period = 97.7 ms. 

0xA783, 

0xA641 

ALM.MAGl =0x8341. 

Alarml i5trueifXACCL_OUT>0.5g. 

0xA93C, 

0xA86E 

ALM_MAG2= 0x3CBE. 

Alarm 2 is true If XACCL OUT < -0.5 g. 
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FRONT VIEW 



Figure 76. 24-Lead Module with Connector Interface 
(ML-24-2) 

Dimensions shown in millimeters 



ORDERING GUIDE 


Model 

Temperature Range 

Package Description 

Package Option 

ADIS16405BMLZ' 

ADIS16405/PCBZ' 

-40T to+105“C 

24-Lead Module with Conneaor Interface 

Interface Board 

ML-24-2 


' Z = Rc^S Compliant Part. 


02009 Analog Dovkot, Inc. All lighu roMrvod. Tradomarki and 
rogiitorod tractemarks aro th« proporty of thair latpactiva ownars. 

007907-0-3/09(0) 


□ ANALOG 
DEVICES 
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C. SIMULINK® MODEL 
1. Overall Model 

ITINYSCOPE Simulation Modj i 
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Figure 59. TINYSCOPE Overall Simulink® Model. 
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Figure 60. Orbital Propagator Block. 

3. Dynamics and Kinematics 
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Figure 61. Spacecraft Dynamics and Kinematics Block. 
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4, 


Environmental Effects 
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Figure 62. Environmental Effects Block. 
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5. 


Disturbance Torques 



Torque 

Solar 


Figure 63. Disturbance Torque Block. 
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6 . 


Attitude Sensors 
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Figure 64. Simulation Attitude Sensor Block. 


D. MATLAB® CODE 
1. Main Script 

%% TINYSCOPE Main Script 
"6 

% Author: 

% LCDR J. Allen Blocker 
% Naval Postgraduate School 
"6 
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% LT Jason D. Tuthill 
% Naval Postgraduate School 
"6 

% Thesis Advisor: 

% Dr. Marcello Romano 
% Naval Postgraduate School 
"6 

%% Format 
clear all 
close all 
clc 

global CONST 
R2D = 180/pi; 

D2R = pi/180; 

%% Set Simulation Conditions 

InitialEulerN =[00 0];%deg 
ReferenceEuler =[00 0];%deg 


Kp = .25; 

Kdl = .15; 

Kd2 = .15; 

Kd3 = .15; 

Ki = .15; 

*** Toggle switches turn the labeled functions on (1) or off (0). 


Tgg_toggle 

= 

0; 

Taero_toggle 

= 

0; 

Tsolar_toggle 

= 

0; 

timeOn 

= 

1; 

taOn 

= 

0; 

cboOn 

= 

0; 

qbnOn 

= 

1; 

qbnmOn 

= 

1; 

rOn 

= 

0; 

hOn 

= 

0; 

e32lOn 

= 

1; 

wbnOn 

= 

1; 

tcOn 

= 

0; 

hsOn 

= 

0; 

wbnfOn 

= 

1; 

biasOn 

= 

1; 

biasfOn 

= 

1; 

pdOn 

= 

1; 

pnOn 

= 

1; 

qbnfOn 

= 

1; 

wbnmOn 

= 

1; 

werrOn 

= 

1; 

berrOn 

= 

1; 

qerrOn 

= 

1; 
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%% Set Constants 
CONST.mu 

CONST.mu_moon = 

CONST.mu_sun = 

CONST.Re 

CONST.Rs 

CONST.J2 

CONST.J3 

CONST.J4 

CONST.SolarPress= 
CONST.SOLARSEC = 
CONST.w_earth 
CONST.Cd 
CONST.Cr 
fleet 

CONST.OmegaDot = 
vance for sun-synch 


398.6004418el2 ; %m''3/s''2 
4.902802953597el2; %m''3/s''2 
1.32 7122E2 0; %m^3/s^2 
6.378137E6; %m 
1.4959787ell; %m 
1.08262668355E-3; % 

-2.53265648533E-6; % 

-1.61962159137E-6; % 

4.51e-6; %N/m^2 
806.81112382429; %TU 
[0;0; .000072 9211585530]; %r/s 
2.5;% 

. 6 ; % 


1.991e-7; %rad/s 


earth radius 
solar radius 
J2 term 
J3 term 
J4 term 

solar wind pressure 

earth rotation 
Coefficient of Drag 
Coefficient of Re¬ 
ascending node ad- 


%% Set Orbital Elements 

%Kep elements meters and radians (a,e,1,W,w,n) 


h_p = 500e3;%m 

h_a = 500e3;%m 


altitude at perigee 
altitude at apogee 


RAAN = 0;%rad 

w = 0;%rad 

TAo = 0;%rad 

Rp = CONST.Re+h_p; %m 

Ra = CONST.Re+h_a; %m 

e = (Ra-Rp) / (Ra+Rp) ; % (m/m) 

a = (Ra+Rp)/2;%m 

ho = sqrt(a*CONST.mu*(l-e^2)); %m“2/s 
mentum 


Right Ascention 
argument of perigee 
true anomaly 
radius of perigee 
radius of apogee 
eccentricity 
semi-major axis 
initial angular mo- 


P = 2*pi*sqrt(a^3/CONST.mu); %sec Orbit Period 

i_sunsynch = acosd ( (CONST . OmegaDot* (1-0*2 )''2*a''(7/2 )).. . 

/(-3/2*sqrt(CONST.mu)*CONST.J2*CONST.Re*2));%eqn 4.47 from Cur¬ 


tis 


i = i_sunsynch*D2R; %deg (rad) 


orbit inclination 


[Ro,Vo] = sv_from_coe(CONST.mu, [ho e RAAN i w TAo]);% initial 
orbital state vector 


%% Set ICs 

ON = DCM(1,-90)*DCM(3,TAo+90)*DCM(1,i*R2D); 
InitialEulerO = DCM2Eul(ON); 

w_BNo = [0;-2*pi/P; 0] ; %rad initial body rates 
w_ON = [0;-2*pi/P;0]; %rad 

% Sensor parameters 
% Gyro 

GYRO_Bias = (3*randn(3,1))*pi/180; % + 3 deg(rad)/sec 

N_ARW = (0.029)*pi/180; 

K_RRW = (0.0002)*pi/180; 
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ARW = N_ARW^2; % angular white noise Variance 

RRW = K_RRW''2/3; % bias variance 

Gg = eye(3) .* (-0.01 + 0.02*rand(3)) +... 

(ones(3,3)-eye(3)).*(-0.0006+0.0012*rand(3)); %percent 

% Magnetometer 
sigMag = 1.25e-7; 

%Mag_bias = (4*randn(3,1))*le-7; % Tesla, +/- 4 mguass 

Gm = eye(3).*(-0.02+0.04*rand(3)) +... 

(ones(3,3)-eye (3)) .*(-0.0028 + 0.0056*rand(3)); %percent 


% Sun Sensor 

51 = [0 20 0]'*pi/180; 

52 = [0 160 0]'*pi/180; 

%Gss = eye(2).*(-0.02+0.04*rand(2)) +... 

% (ones(2,2)-eye(3)).*(-0.0028+0.0056*rand(2)); %percent 
FOV = 0.7; 
sigSS = 0.1; 

J = Bessel(sigSS/3,FOV).*pi/180; 

% Star Tracker 

sigST = 70 /3 /60 /60*pi/180; %arcsec to rad (3sig) 

%Gst = eye(3) .* (-0.02 + 0.04*rand(3)) +... 

% (ones(3,3)-eye(3)).*(-0.0028+0.0056*rand(3)); %percent 

% Kalman Filter 
dt = 0.05; 
t_ekf = dt; 
sig(l) = sqrt(ARW); 
sig(2) = sqrt(RRW); 
sig(3) = sigST; 
sig(4) = sigSS*pi/180; 
sig(5) = sigMag; 

ReferenceOmega = w_ON; 

[qBNo] = Euler_to_Quaternion(InitialEulerN); 

[qBOo] = Euler_to_Quaternion(InitialEulerO) ; 

[ReferenceQuaternion] = Euler_to_Quaternion(ReferenceEuler); 


%sec (20 Hz) model speed 
%sec (20 Hz) ekf speed 
%rad/sec^(1/2), ARW 
%rad/sec^(3/2) , RRW 
%rad. Star Tracker Error 
%rad. Sun Sensor Error 
%tesla, magnetometer error 


%% Run Simulation 

[Spacecraft]= Calculate_6U_Spacecraft; 
J_Matrix = Spacecraft.MOT; 
[density_table] = GetDensity; 

RunTime = P;%sec 
tic 

sim( 'TinyscopeMainModel' ,RunTime); 
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Total_Model_time = toe 

factor = RunTime/Total_Model_time 

DisturbanceTorques.Tgg = Tgg; 

DisturbanceTorques.Taero = Taero; 

DisturbanceTorques.Tsolar = Tsolar; 

SensorMeasurements.ST = q_BNm; 

SensorMeasurements.Gyro = w_BNm; 

SensorMeasurements.bias = bias; 

SensorMeasurements.SSI = ssl; 

SensorMeasurements.SS2 = ss2; 

SensorMeasurements.Mag = squeeze(Bm) 

EKFerror.bias = squeeze(b_err); 

EKFerror.rate = squeeze(w_err); 

EKFerror.quat = q_err; 

EKFerror.cov = Pdiag; 

EKFerror.Pnorm = Pnorm; 

FilterEst.Q = squeeze(q_BNf) 

FilterEst.Gyro = squeeze(w_BNf) 

FilterEst.bias = squeeze(bias_f)'; 

%% Output Results 

% PlotOrbit(R); 

"6 

% PlotEKFerrors(EKFerror,SimTime); 

"6 

% PlotMeasurements(SensorMeasurements,SimTime); 

% PlotFilter(FilterEst,SimTime); 

% PlotModelStats(w_BNs,Bs,SimTime); 

% PlotQuaternion (squeeze(q_BNf)',q_BN,SimTime); 

% PlotDisturbanceTorques (DisturbanceTorques,SimTime); 

% PlotControlTorques (Tcontrol,SimTime); 

% PlotTorques (DisturbanceTorques,Tcontrol,SimTime); 

% PlotEulerAngles (Euler321,SimTime); 

"6 

% PlotQuatError (q_err, SimTime); 
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2 . 


Attitude Matrix 


%#eml 

function att = ATT( quat ) 

%% Making Attitude Matrix 

att = transpose(XT(quat)) * PSI(quat); 

return 


3. XI 


%#eml 

function xi = XT( quat ) 

%% Making Xi Matrix 

xi = [ quat(4)*eye(3) + SKEW(quat(1:3,1)) ; - quat(1:3,1)' ]; 

return 


4. PSI 


%#eml 

function psi = PSI( quat ) 

%% Making Psi Matrix 

psi = [ quat(4)*eye(3) - SKEW(quat(1:3,1)) ; - quat(1:3,1)' ]; 

return 

"6 

5. Skew 


%#eml 

function sk = SKEW( vec ) 

% Check it is a 3x1 or 1x3 vector 

if ~(( (size(vec,l) == 3) && (size(vec,2) == 1) ) || ( (size(vec,l) == 

1) && (size(vec,2) == 3) )) 

disp('not a vector'); 
return 

end 

sk = [0 -vec (3) vec (2); 

vec (3) 0 -vec(1); 

-vec (2) vec(1) 0 ] ; 

return 


6. Quaternion to Euler 

function euler = Quat2Euler( quat ) 
euler=zeros (3,1) ; 

a = 2* (quat (4) *quat (2) -quat (3) *quat (1) ) ; 
if a > 1 
a = 1; 

elseif a < -1 
a = -1; 

end 
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euler(1)=atan2 ( 2*(quat(1)*quat(4)+quat(2)*quat (3)), 1- 
2 * (quat (1) *quat (1) +quat (2) *quat (2) ) ) ; 
euler(2)=asin(a); 

euler(3)=atan2 ( 2*(quat(3)*quat(4)+quat(1)*quat (2)), 1- 
2* (quat (2) *quat (2) +quat (3) *quat (3) ) ) ; 
return 


7. Bessel 

function [J] = Bessel(sig_SS,FOV) 

%% This function creates the spacial noise of the Sun Sensor 

r = linspace(eps,FOV); 
theta = linspace(0,2*pi) 
n = length (r); 

J = zeros(n,n); 

for i = 1:n 

J(i,:) = besseli(rand(n,1),r*r(i))*cos(theta)*sig_SS; 

end 

J(:,l) = J(:,1)./norm(J(:,1)); 

J(l,2:n) = J(1,2:n)./norm(J(1,2:n)); 

return 

"6 


8. Multiplicative Quaternion Extended Kalman Filter 

function [wkl,qkl,biaskl,Pkl] = EKF(wklt,ykl,bs,S,Bkl,B,dt,sig,mflag) 

sig_v = sig(1); 
sig_u = sig (2); 
sig_st = sig(3); 
sig_ss = sig(4); 
sig_mag = sig(5); 

persistent qk biask wk Pk; 

% Initialize States and Measurement 
if isempty(qk) 

qk = [0 0 0 1 ] ' ; 
biask = zeros(3,1); 
wk = wklt; 

sig_n=sqrt(sig_st^2+sig_ss^2+sig_mag^2); 

Pk = 

dt^(1/4)*sig_n^(1/2)*(sig_v^2+2*sig_u*sig_v*dt^(1/2))^(1/4)*eye(6); 

wkl=wk; 
qkl=qk; 
biaskl=biask; 

Pkl=Pk; 
return ; 

end 

MaxST = 1; 
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MaxSS = 3; 
MaxMag = 4; 


Index of last sun sensor 
Index of last magnetometer 


"5 

g, 

0 


SSangles = [20;160]; % Sun Sensors fram angles 

SSaxis = [2; 2]; % Sun Sensor rotation axis l=x, 2=y, 3=z 

%% Propagation 
biaskl = biask; 


Skew_w = SKEW(wk); 

Mag_w = norm(wk); 

psik = (sin(l/2*Mag_w*dt)/Mag_w)*wk; 

Omega = [cos(l/2*Mag_w*dt)*eye(3)-SKEW(psik) psik; 

-psik' cos(l/2*Mag_w*dt) ]; 

qkl = Omega*qk; 

Phi_ll = eye(3)-Skew_w*sin(Mag_w*dt)/Mag_w + Skew_w^2*(l- 
cos(Mag_w*dt))/Mag_w^2; % 7.59b 

Phi_12 = Skew_w*(1-cos(Mag_w*dt))/Mag_w^2 - eye(3)*dt -... 
7 . 59c 

Skew_w^2*(Mag_w*dt-sin(Mag_w*dt))/Mag_w^3; 


Phi_21 

= zeros(3); 

"6 

7 . 59d 



Phi_22 

= eye (3); 

0, 

0 

7 . 59e 



Phi = 

7 . 59a 

[Phi_ll Phi_12; Phi_21 Phi_22]; 

"6 

Gk = [ 

-eye (3) zeros (3); zeros (3) eye(3)]; 


Qk = [ 

(sig_v'2*dt + l/3*sig_u''2*dt''3) *eye (3) 

- (l/2*sig_u''2*dt^2) *eye (3) 

r 

] ; 

- (l/2*sig_u''2*dt''2) *eye (3) 

(sig_u''2*dt) *eye (3) 

Pkl = 

Phi*Pk*Phi'+Gk*Qk*Gk' ; 



%% Update Loop - 

if (sum(mflag) >= 1) 

Att = ATT(qkl); 
delX = zeros (6,1); 
for i = l:MaxMag 

% Compute H matrix for Star Tracker Measurement 
if ( (mflag(i) == 1) && (i <= MaxST) ) 

Xi = XI(qkl); 

H= [ l/2*Xi (1:3, : ) zeros(3,3) ]; 

R = sig_st^2*eye(3); 
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% Gain 

K = Pkl*H'/(H*Pkl*H' + R); 

% Update 

Pkl = (eye(6) - K*H)*Pkl; 

res = ykl(l:3,l) - qkl(l:3,l); 
delX = delX + K*(res-H*delX); 

% Update for Sun Sensor Measurement - 

elseif ( (mflag(i) == 1) && (i <= MaxSS) ) % to max number of 

Sun Sensors 

H = [ DCM(SSaxis(i-MaxST),SSangles(i-MaxST))*SKEW(Att*S) 
zeros(3,3) ] ; 

R = sig_ss^2*eye(3); 

% Gain 

K = (Pkl*H')/(H*Pkl*H' + R); 

% Update 

Pkl = (eye(6) - K*H)*Pkl; 

res = bs(:,i-MaxST) - DCM(SSaxis(i-MaxST),SSangles(i- 
MaxST))*Att*S; 

delX = delX + K*(res-H*delX); 


% Update for Magnetometer Measurement - 

elseif ( (mflag(i) == 1) && (1 <= MaxMag) ) % to max number of 

Magnetometers 

H = [SKEW(Att*B) zeros (3,3)]; 

R = sig_mag''2*eye (3) ; 

% Gain 

K = (Pkl*H')/(H*Pkl*H' + R); 

% Update 

Pkl = (eye(6) - K*H)*Pkl; 

res = Bkl - Att*B; 

delX = delX + K*(res-H*delX); 

end 

end 

qkl = qkl+l/2*XI(qkl)*delX(1:3,:); 
qkl = qnormalize(qkl'*qkl, qkl) ; 
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biaskl = biaskl + delX(4:6,:); 

else 

end 

wkl = wklt - biaskl; 

% Save previous values 
qk = qkl; 
biask = biaskl; 
wk = wkl; 

Pk = Pkl; 

return 

"6 


%% Normalizing routine for quaternions 
function qkl = qnormalize(qnorm,qkl) 
while (qnorm) > 1 

if qnorm > 1 + le-12 

qkl = ((3 + qnorm)/(I + 3*qnorm))*qkl; 
% rescale quaternion to (err^3)/32 

else 

qkl = qkl/sqrt(qnorm); 

% renormalize quaternion 

end 

qnorm = qkl'*qkl; 

end 

return 


%% Make Sun Sensor DCM 
function A = DCM(axis,a) 

A = zeros(3,3); 
switch axis 
case 1 

A = [1 0 0; 0 cosd(a) sind(a); 0 -sind(a) cosd(a)]; 
case 2 

A = [cosd(a) 0 -sind(a); 010; sind(a) 0 cosd(a)]; 
case 3 

A = [cosd(a) sind(a) 0; -sind(a) cosd(a) 0; 001]; 

end 

return 
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